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1. SUMMARY 


In this report, we have described a method for the computation of viscous 
transonic flows over 3-D wings using a zonal approach to treat the viscid- 
inviscid interactions. The chord Reynolds number was considered large and the 
boundary layer was assumed to be predominantly turbulent. For the inviscid 
flow computation, a parabolic coordinate mapping was used in conjunction with 
a finite volume formulation of the conservative full potential equation. A 
new numerical -AFZ scheme was developed for the 3-D Inviscid flow solution to 
replace the SLOR scheme. A special far field asymptotic boundary condition 
was derived that gave more accurate results and better convergence 
performance. In addition, a second order artificial viscosity is used in the 
supersonic zone rendering the computation formally second-order everywhere 
except at the captured shocks. For the 3-D boundary layer calculation, the 
integral method of Ityri ng-Sml th-Stock was extensively modified and made 
suitable for our interaction calculation. The wing thickness effect was taken 
Into account and the viscous wake solution was computed beyond the trailing 
edges. The Interaction calculation was formulated with a set of coupling 
conditions that Included the proper source flux distribution due to the 
surface boundary layers on the wing, the flux jump distribution due to the 
viscous wake and the viscous wake curvature effect. Transpiration boundary 
conditions were used for the inviscid flow boundary conditions for the coupled 
calculations. In addition, a method was devised so that the results of the 
trailing edge strong interaction solution in our 2-D viscous airfoil analysis 
could be adapted for the normal pressure correction near the trailing edge 
region. The wake surface (that is, a fictitious surface imbedded in the 
viscous wake) was floated such that the converged solution coincides with the 
inviscid flow stream surface. A computer program was written to perform 
calculations Including all the above mentioned effects in a fully automated 
manner. A standard case of calculation Is represented by data Inputs of wing 
geometry, section ordinates, freestream Mach number, angle of attack and mean 
chord Reynolds number. Three versions of the source language code of the 
computer program were prepared. One scalar version Is to be used on the IBM- 
3081 computer, and two vectorized versions are to be used for the computers 
Cyber-205 and Cray-lS. For engineering requirements, a typical case of 
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calculation for a viscous solution usually takes about 5 minutes of CPU 
computing time on a Cray-lS computer. Details of the program can be found in 
Volume 2 of this report - "GRUMWING User's Manual". The viscous program has 
been test run on the Lockheed Wing A - a transonic supercritical transport 
wing, and on the Lockheed Wing B - a high subsonic cruise fighter wing. In 
addition, many check calculations have been carried out with the basic AFZ 
inviscid code. These include computations over an 0NERA-M6 wing modified to 
include a large degree of wing twist, and calculations for a more complex 
swept wing designed at Grumman Aerospace Corporation. 
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2. INTRODUCTION 


There has been steady progress In the computation of inviscid transonic 
flows in the decade following the original contribution of Munnan and Cole in 
1971. Recent developments (reviewed in Ref. 1-3) have led to fast, reliable 
and accurate methods for solving the full potential flow equation using ADI 
and multi grid techniques for transonic flows over airfoils, wings and 
relatively complex wing-body combinations. Major advances have also been made 
in the solution of the full Euler equations using multi grid and Runge-Kutta 
techniques. The availability of these codes coupled with advances In 
computers has revolutionized the aerodynamic design process so that 
computational methods currently have a central role in aerodynamic design. 

It is of course well known that viscous effects are important In most 
aerodynamic flows of Interest, particularly at transonic speeds, and must be 
Included In the theoretical formulation to be useful for practical design 
applications. Recent applications of Inviscid codes to transonic wing design 
reported In the papers of Ref. 4 showed the Inadequacies of purely inviscid 
methods for practical wing design at transonic speeds. At the high Reynolds 
numbers of interest In most aerodynamic problems the flow field Is essentially 
Inviscid with viscous effects confined to thin layers near the surface and in 
the wake. In these problems viscous effects can be taken into account through 
a zonal approach based on an Iterative solution of the combined inviscid and 
boundary layer equations, including other local regions of strong vlscld- 
inviscid iteration where necessary. An alternative approach Is that of 
directly solving the full or parabolized approximation of the Reynolds 
averaged equations numerically. While offering the prospect of providing the 
most accurate and general method for solving viscous flow problems, Navler- 
Stokes solvers have not yet led to a practical methods for aerodynamic 
design. This is due to the Navier-Stokes solver's large computing 
requirements and expense and the Inadequacies of the turbulence models 
available for finite difference formulations. In contrast, zonal type methods 
have proved practical for aerodynamic design, particularly when simple 
Integral methods are used for the boundary layer solution. Integral methods 
are attractive for viscld-inviscld interaction calculations not only because 
they greatly reduce the computing requirements but also because they are much 
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more flexible than finite difference methods in adapting turbulence models to 
fit the physics in various subregions of complex turbulent flows (see 
discussion by Kline in Ref. 5). Such methods have proved particularly 
effective for computing viscous transonic flow over airfoils and are widely 
used throughout the aeronautical community. 

The objective of the present work is to develop a zonal type method for 
viscous transonic flow over 3-D wings. Our approach is to apply the same 
viscous/invlscid interaction techniques previously developed for our GRUMFOIL 
code (Ref. 6,7) for airfoils. 

In a zonal type approach the flow is divided into an outer inviscid 
region and thin, inner viscous regions near the airfoil and wake. The viscous 
layers are further subdivided into weak interaction zones comprising most of 
the boundary layer and wake where the standard boundary layer equations apply 
and small local strong interaction regions near shock waves and the trailing 
edge where the boundary layer approximations breaks down due to the presence 
of relatively large normal pressure gradients. The viscous effects on the 
outer Inviscid flow are accounted for through viscous coupling conditions that 
appear as boundary conditions on the Inviscid flow. The ultimate basis for a 
zonal type approach is a large Reynolds number asymptotic limit. To be 
complete in the sense of Including all the leading order terms the matching 
conditions should account for 1) displacement effects on the airfoil or wing, 
2) displacement effects due to the wake and 3) wake curvature effects due to 
the momentum defect across the wake and 4) normal pressure gradient effects in 
the strong interaction zones at trailing edges and shock waves. 

Nearly complete zonal type methods have recently been developed for 
subsonic and transonic flow over airfoils by Melnik (Ref. 6-8) et al , Collyer 
and Lock (Ref. 9), and LeBalleur (Ref. 10). All these methods employed a 
potential flow approximation for the inviscid flow and integral methods for 
the boundary layer solution. They accounted for all the weak interaction 
effects due to the boundary layer and wake but only the GRUMFOIL code of 
Melnik et al accounted for strong Interaction effects at the trailing edge. 
None of the methods accounted for strong interaction effects near shock waves. 

In turbulent flow, the shock wave penetrates into the boundary layer and 
generates large normal pressure gradients leading to a breakdown of the 
standard boundary layer formulation. Although a simple boundary layer 
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description cannot give an accurate prediction of the local details of such an 
Interaction, comparison with experiments (Ref. 6, 8, and 11) indicates that it 
does seem to give a good prediction of the overall thickening of the boundary 
layer across the shock wave and a good prediction of the boundary layer 
integral properties downstream of the Interaction. Fortunately, this seems to 
be sufficient to achieve good predictions of airfoil section characteristics 
even under high-lift supercritical conditions. On the other hand, because 
normal pressure gradient effects near trailing edges can produce significant 
global results through their Influence on the Kutta condition, they should be 
accounted for, particularly when there is a significant degree of rear loading 
on the airfoil. 

A similar method was developed by Whitfield et al (Ref. 11 and 12) using 
an Euler equation solver for the Inviscid flow. Their method only accounted 
for the displacement effect of the boundary layer and wake and ignored the 
wake curvature and strong interaction effects. Although the use of the Euler 
equations for the Inviscid solution might be an improvement for strong shock 
waves (M L ocal > 1*3). the overall method may be less accurate than the 
potential flow methods discussed above because of the neglect of the wake 
curvature and trailing edge interaction effects In their viscous flow 
formulation. For further discussion of developments In zonal methods for 
airfoils see the recent reviews In Ref. 10 and Ref. 12-15. 

Our overall approach to the wing problem addressed in this report has 
been to extend the viscous flow formulation used in the GRUMFOIL code for 
airfoils to three dimensions. For the Inviscid flow, the full potential 
equation for transonic flow over 3-D wings Is employed. A finite volume 
formulation based on Jameson's FL027 code was used to solve the full poten- 
tial equation in conservation form. A modified "C" type mesh was employed to 
allow the wake streamline to be aligned with the grid thus giving a more 
accurate implementation of wake coupling conditions than used in the 2-D 
GRUMFOIL code. A new approximate factorization (AFZ) scheme was formulated to 
accelerate convergence. We also Introduced a more accurate asymptotic far 
field boundary condition and a fully second order artificial viscosity into 
the program. The Inviscid code Is designated FL047. The new laminar and 
turbulent boundary layer method used in the present work has Its origin in 
Stock's method (Ref. 16) but has been extensively modified for Interactive 


5 


calculations. Stock's method uses a fully 3-D integral boundary layer method 
with a lag entrainment method adopted from the work of Green (Ref. 17). 

The viscous flow method employs a fully 3-D form of the wing surface 
boundary layer and wake displacement thickness coupling conditions using a 
transpiration velocity formulation. The wake curvature condition is 
implemented using a quasi -2-D strip formulation that accounts only for the 
pressure variation across the wake that is generated by the streamwise 
component of the curvature of the wake surface. This approximation should be 
useful for the wing flow field under consideration where the primary component 
of the wake curvature is in the streamwise direction. The spanwise component 
of curvature can be expected to be large only near the wing tip, and it is 
only in this region that we showed significant departures from the simple 
strip approximation employed in the present method. We also employ a simple 
2-D strip approximation to justify the use of the same trailing edge 
correction as employed in the GRUMFOIL code. Since gradients in the spanwise 
direction are small compared to the gradients in the streamwise direction, 
this should also be a good approximation except near the wing tip and root 
sections. With these approximations, the present code accounts for all the 
primary viscous interaction aspects except those arising from the shock 
boundary layer interaction. As proved for the airfoil calculations we expect 
this simple boundary layer treatment of shock wave-boundary layer Interactions 
to be adequate for the prediction of wing section characteristics for the 
reasons discussed above. 

Several other methods have recently been developed (Ref. 18-22) for 
computing viscous transonic flow over wings and wing-body combinations. In 
common with the present method, they all employ a modified version of Green's 
lag entrainment method for the boundary layer solution and all, except Samant 
and Wigton (Ref. 22), use a potential flow approximation for the inviscid 
flow. Samant and Wigton (Ref. 22) use the Jameson FL057 (Ref. 23) method for 
solving the full Euler equations in the inviscid flow. Firmin (Ref. 18) and 
Waggoner (Ref. 19) use a transonic small disturbance approach for the inviscid 
flow which is less accurate than the full potential flow approximation 
employed in the present method as well as in Streett's (Ref. 20) and Wigton 
and Yoshihara's (Ref. 21) method. All of the methods, including the present 
one, account for the displacement thickness effect both on the wing surface 
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and in the wake, except for Waggoner's and Samant and Uigton's methods which 
neglect all wake effects. Only Streett's and the present method account for 
the wake curvature effect and only the present method includes the trailing 
edge interaction effect. The corrections to the coupling conditions that 
arise from the local trailing edge interaction solution are particularly 
important for a correct implementation of the wake curvature effect. The 
procedure used for the wake curvature terms in Streett's method, which ignores 
the trailing edge interaction, employs a single mean of the curvatures of the 
upper and lower surfaces of the wake. This undoubtedly leads to a significant 
underestimate of the wake curvature effect. Nevertheless, computations 
presented in Streett's paper Indicated that the effect of the wake curvature 
coupling conditions were fairly large. These results and those obtained 
earlier for airfoils (Ref. 8) clearly demonstrated the importance of including 
the wake curvature terms in the coupling conditions. 

We believe the viscous flow method employed in the present work 
represents the most complete Implementation of a zonal type approach to 
viscous flow over wings achieved to date. 


Acknowledgement. The authors wish to express their appreciation to 
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3. VISCOUS WING THEORY - BACKGROUND 


3.1 INTRODUCTORY REMARKS - INVISCID FLOW 

In the last few years, extensive progress has been made in inviscid 
transonic flow computations. Three dimensional transonic, inviscid flow 
fields can now be computed over wings and wing-fuselage combinations using 
newly developed “finite volume" computer codes which solve the transonic full 
potential equations. For cases where the shocks are not so strong that 
extensive boundary layer separation or entropy production at the shock takes 
place, accurate predictions of lift and drag can be obtained (Ref. 24). Just 
as important as accuracy, these finite volume techniques are flexible, being 
able to handle fairly complex geometries, and can be matched to turbulent 
boundary layer computing methods as reported here and in Ref. 20. 

In Subsection 4.1, we describe a simple but efficient "AFZ" approximate 
factorization scheme which can converge to a solution much faster than 
relaxation methods. The scheme requires no more computer storage than 
relaxation schemes, however, and, as is in those schemes, data is only 
required in a fixed plane- by -plane sequence. Therefore, if necessary, a disk 
or similar mass storage device can efficiently be used. 

Other factorization schemes have been reported for the 3-D transonic 
problem (Ref. 25, 26, and 27). In Ref. 25, two 3-D data arrays must be 
stored, however, and more complicated data manipulations are required due to 
the presence of an additional factor. The "SIP" scheme of Ref. 26 requires 
much more storage and data manipulations. The scheme of Ref. 27 is applied to 
a nonconservative finite difference scheme, rather than a finite volume one. 

In our work, the AFZ scheme is applied to the 3-D conservative, full 
potential, finite volume formulation of Ref. 24. 

Another approach to solving the inviscid equations involves using a 
factorized or a relaxation method for reducing high frequency errors, together 
with a multi grid scheme. Good results, comparable to ours, in terms of 
convergence rates, have been presented in Ref. 28 and 29 for successive line 
over relaxation (SLOR) with multi grid. These schemes require about 30% more 
total storage than ours due to the extra grids used. Also, the use of 
auxiliary storage such as a disk appears to be more difficult due to the 
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required Intergrid transfer operations. However, multlgrld methods may be 
less sensitive to extreme grid stretching than our method. 

An alternative to solving the potential equations for the Invlscld flow 
Involves the Euler equations. These require, however, at least five times as 
much storage and several times as many computations per point per Iteration as 
potential flow methods. Also, although recent progress has been made on 
speeding up convergence (Ref. 30), they require at least several times as many 
Iterations to converge and almost an order of magnitude more total 
computational cost. 

Some of the Initial work on the AFZ method was done with J. Benek and A. 
Jameson (Ref. 31). 

3.2 3-0 BOUNDARY LAYER & VISCOUS MAKE THEORY 

The subject of 3-D boundary layers is a very large and important area of 
research In aerodynamics. It has been of continued interest to aero- 
dynami cists for the past several decades. Examples of recent developments can 
be found In the book “Three-Dimensional Turbulent Boundary Layers" (Ref. 

35). Significant progress has been made on computing attached 3-D boundary 
layers solution on wings. Although 3-D boundary layer separation remains a 
very difficult problem to solve, the phenomenon has become much better 
understood through the study of the development of the characteristics of the 
system of partial differential equations for boundary layers (Ref. 36 and 37). 
The study of boundary layer separation on wings Is particularly Important In 
the transonic flow regime as it Is very often caused by the interaction with 
the invlscld shocks on the upper surface of the wings. 

Methods of solving the 3-D boundary layer equations fall into two main 
categories; finite difference schemes which predict flow properties at each of 
the mesh points across the layer, and integral methods which predict only the 
skin friction and certain Integral properties of the boundary layer. For 
turbulent flow, both categories Involve a good deal of empiricism and both are 
currently underdeveloped particularly when the flow Is separated. There have 
been many Investigations made of 3-D boundary layers on wings. Among those 
using finite difference formulations are the recent studies of Cebeci et al 
(Ref. 38) and McLean et al (Ref. 39). Both studies used eddy viscosity 
modeling. McLean et al made an interaction study by coupling the boundary 
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layer solution to the inviscid flow. The McLean procedure did not account for 
important wake effects and required too much computer time for it to be an 
effective wing design tool. Our 2-D viscous airfoil study indicates that 
there is always a strong wake curvature effect for subsonic trailing edges. 

It Is our belief that this phenomenon Is also important in the 3-D viscous 
wing problem and should be included in the theoretical formulation. 

Integral boundary layer techniques have several advantages over the 
finite difference methods. They are usually easy to apply, the computations 
are fast, and the results are not overly sensitive to the choice of the 
turbulence models used. My ring (Ref. 40) has developed a fairly complete 
momentum Integral prediction method for 3-D boundary layers in incompressible 
flow. The continuity equation (l.e., the entrainment equation) and the two 
Reynolds averaged momentum equations in directions parallel to the surface are 
used to derive a system of first-order partial differential equations for the 
momentum thickness along the mainstream direction, the shape factor, and the 
streamline deviation angle 0 (the angle between the external streamline and 
the limiting streamline at the wall). The global physical quantities (e.g., 
the momentum thickness in the transverse direction) are then implicitly 
dependent upon the deviation angle 0. Myring's work has been extended by 
Smith (Ref. 41) to 3-D compressible turbulent boundary layer flow over an 
insulated surface. To improve the turbulence calculation, Green et al (Ref. 
17) has developed a highly successful prediction method for 2-0 compressible 
turbulent boundary layers and wakes. A one equation lag-entrainment method 
was developed to incorporate nonequilibrium “hi story" effects into the 
integral entrainment method. This method was used in our 2-D GRUMFOIL code. 

The work of Smith and Green et al was combined by Stock (Ref. 16) to 
provide a 3-D extension of the compressible lag entrainment method. Stock's 
method also incorporates an integral method for 3-D laminar boundary layers. 
The laminar boundary layer formulation (Ref. 42) uses five equations, namely, 
continuity, two momentum equations along the surface and their respective 
first-order moment equations, to determine five flow parameters including the 
displacement thickness and the momentum thickness in the mainstream direction. 
The laminar integral is based on one-parameter Falkner-Skan velocity profile 
in the main stream direction and a two-parameter polynomial velocity profile 
for the cross flow. 
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Since Me follow the same formulation as that of Smith and Stock, the 
derivation of the governing equations is not given In this report. However, 
certain modifications and extension of their work which are necessary for 
viscld-inviscid computations are discussed. These Include, the derivation of 
the full metric coefficient taking Into account the wing thickness and the 
variation of the wake shape, the extension of the calculation to 3-0 viscous 
wakes and the derivation of the source distribution from the boundary layer 
solution for the transpiration boundary conditions. 

It Is Important to mention that the resulting 2-D (coordinates along the 
surface) partial differential equations for turbulent flow governing the 
surface shape factor (R), the momentum thickness along the main stream 
direction ( i ) • and the streamline deviation angle e are of hyperbolic type 
with distinct characteristic curves (see Ref. 40 and 41). Similar conclusion 
can also be drawn for the equations for the laminar boundary layer. 

Therefore, the problem for the solution for the system of partial differential 
equation as an Initial value problem is well posed provided proper Initial 
conditions can be prescribed. We shall discuss this more In Subsection 4.2 of 
this report. 

3.3 "ZONAL" APPROACH TO SOLUTIONS OF 3-D VISCOUS FLOW ON WINGS 

From the discussions made in the two previous sections, our Invlscld AFZ 
method for the 3-D potential flow solution and the Integral methods for the 
boundary layer analysis comprise the necessary elements for the vlscid- 
Inviscld interactive analysis for the viscous flow about wings. Our boundary 
layer method follows Stock's (Ref. 16) treatment of the laminar and turbulent 
boundary layer computation on wings with our Improvement that the wing 
thickness effect Is taken Into account. In addition, the complete 3-D viscous 
wake Is computed to take into account the wake thickness effect for viscous 
Interaction. In the present work, the transpiration coupling conditions are 
used for the displacement effect and the wake curvature condition along the 
wake streamline Is determined from the iterative solution of the Inviscid 
flow. In addition, a local 2-D trailing edge strong interaction solution Is 
Included, and Imbedded In the global solution. The derivations of these 
conditions are given in Subsection 4.3 of this report. 


II 


4. NUMERICAL METHODS 


4.1 INVISCID COMPUTATIONAL PROCEDURE 

4.1.1 Basic Equations for 3-D Inviscid Flow 

The basic equation describing irrotational , inviscid compressible flow 
can be written 

a x ( pu) + 3y( pv) + a z (pw) = 0 (4.1.1) 

where 3 X , 3 y , and 3 Z represent partial differentiation with respect to 
Cartesian x, y, and z. It should be pointed out here that the notations used 
in describing the numerical techniques employed in computing the inviscid flow 
(Subsection 4.1) are different from those used for the viscous flow 
formulation (Subsections 4.2 and 4.3). 

When the velocity is related to a potential by 



and the isentropic relation is used for the density 

1 

p= [1 + xM 2 . ,l *q 2 )]^ T • (4.1.3) 

we have a single partial differential equation to solve for the single 
variable, <|>. In the above y Is the ratio of specific heats, q 2 and p are 
normalized by their free stream values and M*, is the Mach number of the free 
stream. 

Equations (4.1.1) to (4.1.3) are transformed to a computational space 
with coordinates (X,Y,Z) by application of coordinate mappings according to 
the method described in Ref. 24. This consists of introducing parabolic 
coordinates (X,Y) In spanwise planes (z = constant) by the transformation 

(X + 1 Y) = /tx-x Q (z)] + i [y-y Q (z)] 

Z- z 
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where x Q (z), y Q (z) define a singular line just Inside the wing leading edge. 
The effect of the transformation Is to unwrap the wing about the singular line 
to form a shallow bump, Y = S()?,Z). The bump Is then removed by a shearing 
transformation, 

X' = 5T, V - r - s(5T,r), V = z 

(see Fig. 1). A far field stretching is then used to define final 
computational variables X, Y, and Z. 

If the transformations are represented by a Jacobian matrix 


V 

3yX 

¥\ 


V 

a y y 

3 z y 

(4.1.4) 

V 

dyZ 

h z J 



Eq (4.1.1) m^y be written In terms of transformed variables: 

L 0 U) a 3 x (phU) + 3y ( phV) + a z (phW) = 0 (4.1.5) 

where 

h = det(H), 

U, V, and W are contravarlant velocities: 



(4.1.6) 


and the physical velocities are given In terms of derivatives of the potential 
In the transformed space: 



(4.1.7) 
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4.1.2 Boundary Conditions 


Initially, freestream conditions were set on the outer boundaries in the 
form <fr = 0 on Y = Y M and = 0 on the surfaces X = X M and X = - X M (Fig. 

1). These surfaces are far behind the wing in the physical space and are 
separated by the wake lY * 0, Z < Zj) shed from the trailing edge of the 
wing. The condition = 0 requires that the streamwise tangential derivative 
of $ be zero at the wake in the far field. All boundaries are at a finite 
distance from the wing in the procedure used here so that it was felt that 
more accurate boundary conditions could be derived by assuming that the 
Prandtl-Glauert equation is valid at these boundaries and solving for the 
reduced velocity potential. The solution can be written formally as; 

* ~ ♦ LIFT + 0 (“j) + 0 


where 


♦LIFT 


_ y 
■ va 


1 1 
w.s. 


a*U,c) 

H 

(z-t) 2 +y 2 


(x - O 


]<U d? 


R » {(x - 5 ) 2 + 3 2 C(z - c) 2 + (y - n) 2 


1/2 

]} 



5,n,c = surface coordinates of the wing (see Fig. 2). 


where the integration is carried out over the wing surface (w.s.) The wing 

thickness integral of the order 0 , and the volume integral of the 

1 R 

order 0 (-y) are neglected in the far field compared to On the 

R p 

boundaries considered here (see Fig. 2), it is seen that, when (y - n) « 1, 

then (x - C) 2 » 0. Also when x “0(5), it follows that on the boundary, that 

lyl >> Ini, irrespective of the values of z and c. A very good approximation 

for — — is then 
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(x - £) 


(X - 0 

— K — 


{(x-e) 2 + e 2 [y 2 + (z - c) 2 


177 * 

]} 


Further expansion to O(^) allows us to express tufj as: 
%ar Field = *LIFT * *1 + *11 


with 


y +b/2 (1 + w^rUJdc 

■ 71 -In u - c ) 2 * y 2 


(4.1.8a) 


and, 


4 -* T 

♦il " Ti J 


dc 


K 9 Ij{Jl)MP{X^)f r 

-b/2 (z-c) +y 


- r(?) } (4.1.8b) 


where 


TT, = 


x{x 2 + e 2 [y 2 + (z - c) 2 ]} 


- 1/2 


"2 
r(c) 




{x 2 { x 2 + e 2 [y 2 - (z-c) 2 ] } 


-3/2 


{ X 2 ♦ eV ♦ (z-t) 2 ] } ' l/2 | 


= U,c)d£ » A<>U) (KUTTA) 


r(c) = # <frU,c)d 5 


It can be shown that the expression for <|>j approaches the relation given by 
K1 anker (Ref. 32) only for x + + « and differs from it otherwise. When x ■*> - » 
0 which Is the correct behavior for the reduced potential In this 
limit. Although the code has been written to Include both the <j>j and <j>jj 
terms in the far field boundary conditions we have so far made extensive runs 
with only the first term included. The results obtained to date Indicate that 
for a value of BOUND ■ .9 (i.e., on the boundary about 4 chords from the 
wing), using <j>j alone Is quite sufficient. The computed lift results are 
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consistently higher than those satisfying free stream velocity at the 
boundary, indicating the importance of the effect of the far field boundary 
conditions. Also, the convergence rate is better than that achieved using the 
original freestream conditions previously discussed. Details of computed 
results will be discussed in Subsection 5.1 of this report. 

The boundary condition on the surface Z = Z M (Fig. 1) is <}> = 0. A 
symmetry condition is used on Z * 0 (i.e., <|> Is reflected) since the flow is 
symmetric there. 

The surface Y = 0 in the computational space contains the wing for -Xy >E< 
< X < X T £ and Z < Zy. The velocity tangency condition V = 0 is Imposed here 
by setting = 0 on the wing surface. The wing trailing edge is mapped to 
two lines, X = Xy E and X = -Xy E , Z < Zy in this plane, so that the wake from 

the wing trailing edge Is mapped to two parts of the surface: X > x TE and X < 

- Xj E , Z < Zy. Across these surfaces the normal velocity V and the pressure 
(p) are forced to be continuous so that both p and V must match at the two 
computational points corresponding to the same physical location. Since p is 
a function of q 2 alone in potential flow, and V and W are small compared to U 
(I.e., the coordinate surface Is approximately aligned with the free stream) 
continuity of pressure is, to a good approximation, enforced by requiring 
continuity of U. This, in turn, requires continuity of fy. Since, for a 

lifting wing <f> is not continuous across the trailing edge it is not continuous 

across these surfaces. 

There are two options for positioning the wake in the computer code. The 
standard one Involves keeping It fixed as the solution converges so that, in 
general, there will be (continuous) flow through It. The second option 
Involves allowing the wake to follow the local velocity. With the latter the 
wake Is still mapped to part of the Y = 0 surface and the mesh is Iteratively 
adapted to the flow. This mapping requires that the height (Y value) of the 
wake be a single valued function of X and Z, which precludes its rolling up. 
The equations for this option are described in Subsection 4.1.5. Off the wing 
tip in the Y = 0 plane (Z > Zy, all X) continuity of V and is enforced but 
no discontinuity of $ Is allowed. The potential Is then matched at points 
which coincide in the physical plane. 

Two exceptional lines occur In the Y = 0 plane where special conditions 
are applied. The first is the line corresponding to the mapping singularity X 
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2 2 

* 0 for Zj < Z < Z M . There, 3^$ + 3y 4> * 0 is enforced. The other 
corresponds to the wing tip, Z ■ Z T , -Xj E < X < Xj E . There the physical 
vortex sheet would start to roll up. Since the computation constrains the 
sheet to start behind the wing and disallows roll up, an additional condition 
is required to avoid a singularity in V. It should be pointed out that this 
flat vortex sheet constraint is not inherent in the assumption of potential 
flow. The existence of a potential is consistent with the appearance of 
rolled up trailing vortex sheets or vortices at wing tips. Both these 
phenomenon can be modelled in full potential flow. The condition we impose 

(at the wing tip in the code) is that 3y<f> = 0 all along this line, which 

insures that V is a local maximum. Thus far, the use of this condition has 
proved necessary in flows with highly loaded wing tip sections. 

4.1.3 Discretization 

The finite volume scheme of Ref. 24 is used to discretize Eq (4.1.5). 

Two staggered grids are used (see Fig. 3a). On one, p, u, v, w, U, V, W, H 
and h are defined; on the other, x, y, z, X, Y, Z, $, and L( «j») are defined. 

The nodes of each grid are at the centers of the cells of the other, and a box 

scheme is used to compute the derivatives in Eq (4.1.4), (4.1.5), and 

(4.1.7a). For a variable f, 

V ’i+lj+l.k+l = [ f i+l,j+l,k+l + f i+l,j ,k+l + f i+l,j+l + f i+l,j,k 

“ f i ,j+l,k+l “ f i,j,k+l ' f i ,j+l,k * f i,j,kJ /4 

Similar expressions can be written for the Y and Z derivatives. 

The indices i, j, k are used for the computational coordinates X, Y, and 
Z which are normalized such that 


AX ■ AY « AZ = 1 


Using these relations, all the partial derivatives in Eq (4.1.4), (4.1.5) and 
(4.1.7) are replaced by finite differences and the equations become finite 
difference equations. 
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The wing surface bisects the flux balance cell such that four cell 

corners are above the surface and four image corners below. The values of U 

and W below the surface occuring in the discrete version of Eq (4.1.5) are set 
equal to those above, and the values of V set equal to the negative of the 
values above. This is described, in 2-D, in Fig. 3b. Across the wake the 
same flux balance is used as in the interior. The values of U, V, and W at 
the image corners on one side are taken to be the values above the surface on 

the other side. This defines the equations on one side of the wake in the 

computational plane and ensures that V is continuous to second order at 
convergence. On the other side of the wake values of <|> are used equal to 
values at corresponding points on the first side, with a constant 
discontinuity (in each Z plane). The value of this discontinuity is 
determined at the trailing edge points in each plane. This enforces 
continuity of pressure as described above. 

The same procedure is used beyond the wing tip where there is no 
discontinuity. The far field boundary conditions are evaluated by numerically 
integrating Eq (4.1.8) with the jump in potential given at the trailing edge 
of each Z * constant section for Z < Zj. 

The basic finite volume scheme leads to an odd-even decoupling of 
solutions. If a small term is added to the left hand side ( L 0 ( <j>) ) of Eq 
(4.1.5), this problem can be eliminated. First, L(<j>) is expanded: 

2 2 2 

L( <j>) * A 3Jj<|> + B 3y<j> + C 3^ + other terms. 

The coefficients A, B, and C can be found from the nonconservative quasi 11 near 
formulation: 


A = pMgjj - U 2 /a 2 ) 

B * P h (9 2>2 - V 2 /a 2 ) 

C = ph(g 3>3 - W 2 /a 2 ) 

where g^ j are the diagonal elements of (H T H) -1 and a is the local speed of 
sound. The recoupling terms are then 
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L rec U) “ 6 xy (A+B)6 XY * + 6 YZ (B+C)6 YZ <fr + 5 ZX (C+A) ^X**’ 


+ 5 xyz (A+B+C) 6 xyz ^ 


Here, 


6 XY f L.l .J- uJ- = ^ f i+lJ+l,k+l ' f 1 ,j+l,k+l ‘ f i+l,j ,k+l + f 1 ,3 ,k+l 

+ f 1+l,j+l,k " f 1,j+l,k " f i+l,j ,k + f 1,j,k^ /2 

etc, and 

6 XYZ f “ t f 1+l,J+l,k+l ' f 1,j+l,k+l * f 1+l,j,k+l + f i,j,k+l 

^i+l,j+l,k + ^i,j+l,k + ^i+l,j,k " f i,j,k^ 

It can be seen that L rec (<|>) Is of second order with respect to the terras in 
L q ( <j>) . 

For stability in the supersonic zone, either a first order or a second 
order artificial viscosity can be used. The first order form is less accurate 
but results In faster, more reliable convergence to a solution in general. 

For the first order form a switching function is first defined: 

s = max(0, 1 - M 2 /M 2 ) 
c 

where 

H 2 . , V 

and ^ Is a switching Mach number, slightly less than 1 (*> .95). Then, 

n _ #.,2 , 2 , UV , . . WU . .x 

P - y(d ^ $ XY $ + zX^ 

n - fV 2 A X VW X * + UV * *\ 

Q - y(V Sy 4> + -jp 6 Y 2 ( i ) + 7 - 5 xy $> 
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R 3 + -jp- + Tf” 

2 2 

v = cshp/a 

where c Is a constant and the second order difference operators are as defined 
previously, with, in addition, 

6 x f = f i+l,j,k " 2f i,j,k + f i-l,j,k 


etc. 

In terms of these quantities, for flow with U and VI positive and V 
negative, the artificial viscosity term is 

W * P i,j,k ‘ P 1-l,j,k + Q 1,j,k " Q i ,j-l,k + W i ,j ,k “ W i,j,k-1 

(a negative V implies flow in the increasing j direction). 

This is first order compared to L 0 ( <(>) . If this artificial viscosity 
option is chosen, the complete equation to be solved for the potential, <|>, at 
each point is 


+ L rec ( * ) + L v1s (4>) 3 0 

In order to achieve second order accuracy in the supersonic zone L v ^ s «|> is 
replaced by (again for U and W positive and V negative); 

‘-vis* “ (P i,j,k ' e 1,j,k p 1-l,j,k ) " (p 1-l,j ,k “ e i-2J,k P 1-2,j,k ) 

+ similar terms for Q and R. 

The local parameter ^ Is Included for stability at shocks, = 

1 forces the scheme with L v j s <|> added to be second order, while ^ - 0 
forces L v .j s $ to the first order form defined previously. In the work of Ref. 
33 it is shown that the scheme had to be first order near shocks for 
stability. Hence, the form 
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e tJ»k * max(0 * ^ “ v 2 " v 3 ^ ^ 

Is used, where v 2 and v 3 are adjustable parameters and x = ^ k - ^ k 

Is used to detect the shock. This Is small In the smooth parts of the flow 
and large near the shocks. To use the first-order option in the code, v 2 = 1 
and \»2 = 0. For the second-order option, v 2 = 0 and V3 = 4. We have found 
significant differences in the flow field In the supersonic zone between the 
first and second order schemes. For the second order scheme, the shocks are 
captured over fewer mesh points. Results using the first and second order 
schemes will be compared in Subsection 5.1 of this report. 

4.1.4 AFZ Scheme 

The basic Idea behind the AFZ scheme Is to solve a set of Implicit 
equations In each constant Z-plane, for a correction to <t>. Since no Implicit 
equations are solved in the Z direction, the 3-D array of 4> values can be 
stored on disk and only several X-Y planes of data need be stored In the 
computer at any one time. 

During the Iteration of sweep number n+1, when updating plane number k, 
values of t are available corresponding to Iteration n+1 In plane k-1 since It 
has just been updated. Only unupdated (level n) values are available In 
planes k and k+1. It was decided to make use of the available updated values 
In computing the residual at plane k, as is done In success! ve-llne- 
overrelaxatlon (SLOR) where updated values In the previous line are used and 
Implicit equations are solved along each line. It Is most efficient to 
compute contravarlant velocities once each Iteration between each plane and 
use them In computing L 0 (<j>) for both planes on either side. To directly 
Incorporate updated values of <t> Into each L 0 ( 4>) computation, It would be 
necessary to abandon this approach and recompute these velocities, using them 
only once for each L 0 ( 4>) computation. This would almost double the number of 
required calculations. An alternative which was chosen involved adding the 
correction multiplied by an appropriate constant to L 0 (<J>) computed using only 
old values at each plane. Our scheme thus had the form 


Nxy a otuiL ( 4> n ) + auCE z 6<J> n , 


(4.1.9) 


<j> n+ * * ,j, n + 5^, n 
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where N XY is an operator in the XY plane, a may be an operator, w is a 
relaxation factor and E~ is the shifting operator: 


If L( 4> n ) is a Laplacian, 

L l (*) * A<S^4> + B6y<t> + 

then choosing (T = C would make the right hand side, at plane k, equal to 

cm[A«^ + B6 2 $J + C(£j- 2*J + ^ +1 )] 

which Is Guo ( <j>) computed using the latest available values of <j>. For our 
nonlinear case we used a similar approach for determining (T . 

The operator N XY was chosen from successful 2-0 ADI schemes (Ref. 34): 
N XY * “ ®x^X^°Y ” ^ y B 5 y ) 

The values of A and B are given by the expansion of L 0 (<(>), described 
previously. Also, and a Y are numbers at low speed, and, following the 
approach of Ref. 34, become operators for supersonic flow. They also Include 
a part to approximate L v ^ s (<|>) In addition to L 0 ($). 

For flow in the ±X, +Y direction, 

a * a Q + + <*2^ Y 

where 

a o = nax(p 0 ,3 1 ) , 

3 0 = Pj^ Az r(Z) max(0, 1 - q 2 /a 2 ) 

Bj * 2 Pj r(Z) C 

= P 4 P x az r(Z) min (l,q 4 /a 4 ) f (U) 
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og = P 5 P x az r(Z) min (1, q 4 /a 4 ) . 

In the above, Az Is the difference In physical z values between planes, f(U) 
is 1 except near the point where U changes sign, C Is an expansion coefficient 
described previously, P^ - P5 are constants for each iteration and 


r(Z) 


1 + Z ^ Z wing tip * 


We also use 

C ~ = P 2 a o 

except at the root plane (Z = 0), where 

Z * P 2 o 0 /l.l . 

The relaxation factor w Is given by 

a = P 3 (1 + l/r( Z) ) . 

Also, to approximate first order artificial viscosity terms 

.+ .,2,2 
djl 8 a • Ojj 11 U Sjj 

.+ u 2 -2 
oiy s a - iy |i V 0 y 

where u is defined previously. To approximate second order artificial 
viscosity, we use 

0 ^ = a - &+ Ui 2 (l + e)6 Z 
cty * a - 6y yV Z (l - e)6 Z 
where e is defined above. 

The form for Oq was chosen to become small when the local Mach number 
approached 1, which is required for stability. However, modal analysis showed 
that, in far field regions, where the grid is stretched in X and Y, but not in 
Z, there is a lower bound on Oq for stability (for elliptic flow). The 
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functions 3 0 and 3j^ represent these two competing requirements. The 
parameters and were designed to become small for low speed flow so that 
o would then be a number instead of an operator. Finally, it was found that 
the minimum a for stability increased with increasing Z. The function r(Z) 
was accordingly included to improve convergence. The use of l/r(Z) in the 
formula for u> results in a relative under-relaxation for larger values of Z, 
which further increases stability and improves convergence. 

The solution sequence consists of a set of cycles. In each cycle, there 
are K sweeps through the field, in which values of are cycled and set equal 
to P°, eP°, e 2 P°, .... e K-1 P° . All other parameters are kept constant. 

Best results were found for K = 4 and e * 1/2, for fine grid transonic 
calculations. In each sweep, the factorized equations (4.1.9) are solved 
plane-by-plane starting from the wing root, for corrections, 6<j», which are 
added to <)>. In each plane, first, 

( oy - 6^A6^)$ - auiL(<J> n ) + auCE^ti^ (4.1.10) 

is solved row-by-row for a temporary 2-D variable $ . Then, 

(c^ - 6B6y) 6<j> = ^ (4.1.11) 

is solved column-by-column for the correction, 5<j>. The third order operator 
terms in ay are all in one direction, since the Y-velocity always has the same 
sign so that there are four diagonals in the equivalent matrix Eq (4.1.11) for 
5<t>, instead of the three that would occur if only o were used Instead of ay. 
Thus, a four-diagonal solver, which is only slightly more complicated than a 
three-diagonal one is required. For the solution of (4.1.10), a five-diagonal 
solver is required since, in the mapped plane, the flow changes direction 
along the line being solved, and derivatives in both directions occur. Again, 
this does not require much more computation than a tri diagonal solution. 

4.1.5 Fitted Wake Calculation 

In each span station (z-const) a "C" mesh is generated. The inner-most 
coordinate line partly lies on the wing surface. The rest of the line lies 
between the trailing edge and the downstream boundary (see Fig. 4). 
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The grid points are labeled with indices i, j, and k corresponding to 
computational X, Y, and Z coordinates. Our concern will be lines B-A and D- 
E. We will want them to follow the wake streamline from the trailing edge 
(t.e.) to the downstream boundary. Points along these lines will be denoted 
(x,y) = (x| k , yj k ) or (xj k , yT k ) with z = z^, independent of i. We 
define i = 1 at the t.e. and i = N downstream. Points on line B-A are denoted 
by (+) and on line D-E by (-). 

We set xj k s k * and determine 5^ k = y| k - y^ k from a computed 
wake thickness. Then, denoting the mean value 

y i.k ■ W,k + *;,k> - 


and 


x 1,k s x 1,k = x i,k» 

the condition that the median wake line follow the wake streamline is that 


v 1,k = y i+l,k ~ y l,k 
u i,k x 1+l,k " x i,k 


(4.1.12) 


where uj ^ and v^ ^ are mean physical velocities defined in plane k: 




• 7 K,. 


u t,k> 


v i,k *7 lv i,k + v i,k> 

During each Iteration, for which new values of the potential function, $, 
are computed, a displacement Ay-j ^ is computed such that Eq (4.1.12) is 
satisfied. These displacements are then added to each node corresponding to 
the same i and k for the next potential Iteration. The formula used for 
computing the Ay^'s is 

^i+l.k = “ l ^7J (X 1H,k • x 1,k> - ( *1H,k - .k - ^1,k>J 
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where u Is a relaxation factor (presently set to 1). 

The physical velocities u* k v| k are computed from the metric and 
velocity potential 4 values at surrounding points: 

u = " <s x y6 Y < ^ h + cos a 

v = (-6yXS x $ + 6 x x6 y <p)/h + sin a , 


h = fiyyfi^x - 6 x y6yX. 


In the above, for u^ k and v^ k 


++ 


V 8 7 (F ll, k+1 " F i,k 

SyF = T (F + + F ++ - 

Y 1 U i+1, k M,k 


, the differences are computed by 

+ F + - F + 1 
p 1+l,k M,k'» 

F - F + ) 
f i+l,k M,k } 


where superscripts (++) denote quantities on the line above those denoted by 
(+). Similar equations with (+) changed to (-) are used to compute 
u- k , v; k . The velocities are normalized to unity in the free stream and a 
denotes the angle of attack. 

During each iteration a boundary value problem is solved (approximately) 
for the potential <j>, regarding the boundary lines B-A and D-E as fixed. The 
corrections are then added to the coordinate so that these lines follow the 
t.e. streamline, which moves from iteration to Iteration. The boundary 
conditions on <j> are that pressure Is equal on the top and bottom of the lines 
and normal flux through the lines is conserved. These are boundary conditions 
since the lines are mapped to two segments of the outer boundary of the 
computational domain (Fig. 4). Hence, after each <p iteration fluid is flowing 
through the lines. This is then corrected by moving the coordinate system 
using the computed corrections. 

The wake lines bisect a set of cells which have velocities defined on the 
corners, and values of <p and x and y defined in the centers (Fig. 5). To 
second order in cell mesh size, balance of the fluxes in the cells requies 
that the normal flux be continuous through the wake line. This flux balance 
relation is just what is imposed in the interior of the domain. Hence, by 
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extending it to the boundary cells the velocity boundary condition is enforced 
along line A-B. 

The other condition is a force balance, imposed by specifying $ along 
line D-E. This Is used to determine $ s , the derivative along the line. The 
isentropic equation for the pressure, which is used to set a relation for <j» s , 
is 

Y / Y- 1 

p = 1 [1 + 1m* (1-q 2 )] (4.1.13) 

K 

where 

2 2 

q 2 = u 2 + v 2 = U n + qjj’) + (* $ + q*) (4.1.14) 

In the above, q“ are the normal and tangential components of the free 
stream velocity with respect to the line. 

A pressure (pj ^) is first computed along A-B. Then, a desired 
pressure, ( pT . ) is computed for each point along D-E, using a force balance 

1 oo 

relation. Also, a normal velocity, q n = <|> n + q n , Is also computed at each 
point. Equations (4.1.13) and (4.1.14) are then solved for a desired «j>. 

At present, we are only implementing a simplified version of this 
relation: if q 2 « q 2 , then 

Y / Y~1 

p “ "T L l + (1 “ q s )J + 0 (q n /q s ) (4.1.15) 

Also, in this approximation the flow on the upper and lower sides of the wake 
are in the x-y plane. Hence the curvature in the flow direction Is the same 
on both sides and the pressure relation, neglecting viscous effects, is just 

Pl,k = Pi.k 


Using approximation (4.1.15) we then have 


q 


+ 

s 


1 »k 



1,k 
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or 


♦: = * s 

s i ,k s i ,k 

Thus, «j)T k * 4 k + r k ’ wliere r k 1s a constant computed at the t.e. in 
each plane after each iteration. 

Presently, we are implementing a curvature calculation with the 
assumption that the flow is entirely in each x-y plane. Although the inviscid 
pressure balance relation does not require this curvature effect, the viscous 
corrections do. At each point along the wake midway between the coordinate 
points we have the unit normal (see Fig. 6a): 

n i+ V2 »k s { ‘ (y i+l,k‘ y i.k*’ (x i+l,k‘ x i,k ))/DS i-*V 2 ,k 


where 


DS i4/ 2 ,k * ^ x i+l,k ‘ x i ,k ) + (y i+l,k 
The radius of curvature, R, can be defined by 

A$i k 2 

+ 0 <“ > 

where 




^i ,k = ^ x i+l,k ' x i-l,k )2 + (y i+l,k - y 1,k )2 J ^ 

A A 

and A0^ k is the angle between normals n^ k n^j, k : 

A^.k * cos 1 (n^ >k • n ii/ 2 , 1 ^* 

Although we have not implemented a general curvature calculation, which 
would be required when the flow has a large spanwise velocity component, we 
give the formulation here. Considering two lines, 

(i ,k-l) ♦ (i ,k+l) 

and 
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(1-1, k) (1+1, k) 


we can define a vector normal to the wake surface (see Fig. 6b): 

n 1,k = (r i+l,k " r 1,k-l J x (r 1+l,k ‘ r 1-l,k 

This Is second order accurate If the mesh Is smoothly varying. 

We will define curvature In a coordinate system consisting of n^ k , 

6 1 r 1,k 5 (r 1+l,k " r 1 ,k ) * and % ,k 5 n 1,k x 6 1 r 1,k* 

We now take normals, as before, midway between points 1+1, k and 1,k, and 
between i ,k and 1-1, k: 

n i+l/2 »k = (r i+l,k ’ r 1, ) x (r i,k+l ' r 1,k-l ) 

n 1- V 2 »k " (r i,k " r 1-l,k ) x (r 1,k+l " r 1,k-l ) 

The scalar product then gives the relative angle at 1,k for radial lines 
Intersecting ( 1+1 ,k) (i-l,k): 

^i.k s cos 1 [* ”iJL/ 2 ,k ‘ n 1-V 2 ,k )/(|n H l / 2 ,k ,|n 1-V2,k ,) l 

The distance between the mid points, 

^i.k = |r 1+l,k “ r 1-l,k l/2 

then Is used to get the radius of curvature: 



AS 

A8. 


1 


0 (A9 2 ) 
1 ,k 


1 / 

The radius of curvature along the line (1,k+l) ♦ (1,k-l), R. k can then be 
computed in exactly the same way by interchanging the two lines. The radius 
of curvature In the normal direction (S^^) can then be computed: 
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where 




cos e 1jk )/sm 6 ( >k 


e i,k = cos 1 t 6 1 r 1 ,k ' (r 1. 


k+1 


r i,k-l )/(|6 1 


'l,k ll 'i,K + l 


' r 1,k-l 1) 


I 

We have tested the AFZ code with the fitted trailing wake on a number of 
wings. A typical transport type configuration in shown in Fig. 7a. Shown In 
the figure Is a 3-D view of the wing and the Isobars on the upper surface of 
the wing for an angle of attack of 1.5° (upper half of the figure) and 0° 
(lower half). Both cases are at a free stream Mach number of .82. At 1.5° 
angle of attack one can see two shocks at the root section which merge at 
about midspan. Figure 7b shows a 3-D view of the wing and the trailing 
wake. Figure 7c compares the convergence of lift for calculations with and 
without the fitted wake. The figure shows the same convergence rates while 
the final lifts are slightly different. For this case (freestream Mach number 
of .82 and angle of attack of 1.5°) the converged C L * .588 for the fitted 
wake calculation while the calculation without tracking the wake gave C L = 
.590. 


4.2 SOLUTION OF 3-D BOUNDARY LAYER & WAKE 

In this section, we discuss briefly the numerical procedure for obtaining 
the solution of the 3-D boundary layer and wake by the integral method of 
Myrlng-Smlth-Stock. We shall stress the modifications and improvements 
whenever they occur but neglect the detailed derivation of the governing 
equations. Since the solution procedures for the laminar and turbulent 
boundary layers are Identical, only the latter Is discussed. The governing 
differential equations for the Integral laminar boundary layers can be found 
In the report of Stock (Ref. 42). The notations for the coordinates used In 
this section and the subsequent ones In this report are redefined and 
different from those used in Subsection 4.1 of the Invlscid flows. 

4.2.1 Coordinate System & Governing Equations 

A curvilinear nonorthogonal surface coordinate system (x,y) coinciding 
with the wing surface and the wake Is chosen for the differential equations. 
The coordinate x is along the chord direction while the coordinate y is along 
the constant percentage chord line as shown In Fig. 8. Shown also are the 
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velocity components u, v along x, y, respectively, x is the angle between x, 
y axes at any point on the surface, o Is the angle between the external 
streamline and x axis. If the Cartesian coordinates (X,Y,T) of the wing and 
the wake are given, the unique transformation 


X = X(x,y) 

Y = Y(x,y) (4.2.1) 

T * T(x,y) 


exists, from which the metric coefficients, h 1# h 2 , g of the curvilinear 
system can be derived, where, 



Further reduction suitable for numerical Implementation leads to 


h 


1 = 




9 8 



(4.2.3) 


3X 

The term -gy- is related to the sweep of the coordinate line y while jQ. and 

3X 

3j 

•gy are the slopes of the wing/wake surface. Since the boundary layer 
equations will Involve x, y derivatives of h^, h 2 , and g, the curvature terms 
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1' t can affect significantly the solution of the boundary-layer 

3X^ 3Y^ 

equations. 

Denoting the velocity components along the nonorthogonal curvilinear 
coordinate system (x,y,z) (z being normal to the wing/wake surface) as (u,v,w) 
(see Fig. 8 ), respectively, the boundary-layer momentum integral equations as 
given by My ring are, along the x direction. 
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and along the y direction. 


(4.2.4) 
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(4.2.5) 


where M is the Mach number at the edge of the boundary layer and and c f 2 
are the skin-friction coefficients in the x and y directions respectively. 
The velocity components in the x, y directions at the edge of the boundary 
layer are denoted by u^ v^ and the resultant velocity at the boundary-layer 
edge is denoted by u e where 

2 2 2 2 q 

u e = U 1 + V 1 + Kfo Vl * 
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The momentum integral thicknesses 0 31 , 0 12 , and q^ 2 and the mass Integral 
thicknesses ^ and together with the quantities k^ k 2 , k 3 , j^, * 3 and 
q as functions of the metric coefficients h^, h 2 , g are given in Appendix A of 
Ref. 41. 

The entrainment or the integral form of the continuity equation assumes 
the form. 


^q tlx fe” ( u l 6 " Vl^ + |y ( v l 6 " “e^M 

_ 1 r U l 36 V 1 35 w i . F 
u^ ^ 3X + T£ 3y ’ w lJ F 


(4.2.6) 


where p e is the density at the boundary-layer edge and 6 denotes the boundary 
layer thickness. F is the entrainment coefficient, l.e., the nondlmenslonal 
rate of change of mass flow in the boundary layer. The entrainment equation 
is used to provide an independent relation in this procedure by prescribing 
the entrainment coefficient as a function of the local boundary layer 
properties. In the lag entrainment method employed here, this Information is 
supplied through a separate differential equation for the entrainment 
coefficient F derived from an approximation to the turbulent kinetic energy 
equation. 


From the entrainment equation, one can derive a useful expression for the 
source velocity distribution, that is, the normal velocity outflow due to the 
presence of the boundary layer. 


1 r 3 r p e q Vl* a r p e 
m ‘ iyi [ ax t— Rj-J + 
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(4.2.7) 


Another useful form of the entrainment equation is, 

3 , p e qv l*\ .3 r p e q Vl w a r p e qU eV 

_ l _ 1 _ j + -g- (_ j - L— ^ j + -§y l— ^ 


(4.2.8) 


which allows the displacement thickness 6* to be computed when the Integral 
quantities and a 2 are known. 

Equations (4.2.4), (4.2.5), and (4.2.6) form the foundation for the basic 
equations, and In addition, we adopt along the direction of the stream tube 
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lag-entrainment model equation of Green (Ref. 17) for the entrainment 
coefficient F. However, for a closed system, the total number of unknowns 
must be reduced to four. The integral thicknesses e^, e^, Q^, 0£ 2 A 1 » *2 
(i.e., the unknowns) can be expressed in terms of the corresponding 
expressions denoted by 8 ^, e 12 , e 2 i> e 22 5 i» 5 2’ * and ® • These 
quantities, except for X and a , are defined in terms of velocity components 
in the directions along and normal to the local external streamlines. The 
quantity a is the angle between the x-axis and the external streamline and 
X is the angle between the x and y axes. Then by further introducing an 
empirical cross flow velocity profile, we finally can reduce the system of 
four equations for the four unknowns, namely, fl , the equivalent incompress- 
ible shape factor, 8 ^, the momentum thickness defined with the main stream 
direction velocity, F, the entrainment coefficient and y, a cross flow 
boundary layer parameter equal to or implicitly related to the wall limiting 
streamline angle 8 relative to the external streamline, dependent upon whether 
the empirical cross flow profile of Mager (Ref. 43) or Johnston (Ref. 44) is 
used. 

The equations are completed by an expression relating skin friction, Cf, 
to the boundary layer variables 8^, fl , and the external conditions. First, 
we note that the components C fl and C f2 along x and y directions can be 
expressed in terms of the skin friction magnitude, Cf and the angles 
X, a, and 8 by the relations. 


r r f sin(X - a) - cos( X - a) tan ( 8)1 

u fl “ °f l : — ~ J 


sin X 


r r jSin(o) + cos(a) tan(8)i 
« f 1 1 ' 


(4.2.9) 


Following Smith we use the Ludwieg-Tillmann relation for the skin 
friction magnitude, C f , modified for compressible flow according to Eckert's 
reference temperature concept, viz. 


P a u Q 0 n “*268 

C f - .246 (JLe“) 



(4.2.10) 
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where T*/T e = 1 + .13 M 2 for adiabatic flow in air and the coefficient of 
viscosity v* is evaluated at the temperature T* by 



The final form of Eq (4.2.4), (4.2.5), and (4.2.6) together with the lag- 
entrainment differential equation for the entrainment coefficient F, can be 
expressed in the form. 


A au 
A ij 3x 



, 3u 

J 1j W 


« C j 


(4.2.11) 


and u^ = 

The characteristics of the quasi linear system of four partial differential 
equations and their solution are discussed in the following subsections. 

4.2.2 Characteristics & Compatibility Equations 

The equation for the entrainment coefficient F in Eq (4.2.11) is of the 

form, 


U 11 

fl 

Y 


where 
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The x and y derivatives of F at any point can be computed and are decoupled 
from the derivatives of e^, R , and y. Furthermore, the characteristic 
curve for F is the external invlscid streamline. Equations for e^, R , and 
Y are coupled, viz. 
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, a^j * b^j * ^ *0 for i = 1,2,3. Multiplying each 


with 




component of Eq (4.2-12) by and summing, we have, 

( y l a H + u 2 a 21 + y 3 a 31 ^ ^"ax + X 3y ) 

+ ( ^i a i2 + u 2 a 22 + M 3 a 32^3l + X 17^ (4.2.13) 

+ ( + u^ + ^3833) OgJ + x I7) s S*1 + y 2^ + 
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X = 


w l b ll * ^21 * ^3^31 = y l b 12 * y 2 b 22 * y 3 b 32 = * 1*12 * 
y l a ll + w 2 a 21 + y 3 b 31 y l a 12 + w 2 a 22 + y 3 a 32 y l a 13 + 


y 2 b 23 + ^33 
y 2 a 23 + y 3 a 33 


(4.2.14) 


The system of homogeneous equations Eq (4.2.14) has nontrivial solution for 
Up U 2 » a nd U 3 if and only if, 


Detlb^ - a^x| = 0 (4.2.15) 

The solution for X from the cubic equation has been studied rather extensively 
(Ref. 40) In conjunction with the development of the integral methods of 3-D 
turbulent boundary layers. At any point (x,y), X is the tangent of the angle 
between the characteristic line and the x-axis. For incompressible flow, the 
form of the characteristic equation is particularly simple as shown by 
Myring. It was shown that for a well behaved boundary layer Eq (4.2.15) 
possesses three distinct real roots, therefore the system of differential 
equations Is totally hyperpolic. The three characteristic lines lie between 
the external streamline and the limiting wall streamline of the boundary layer 
flow. Numerical solutions of Eq (4.2.15) carried out by us have shown that 
this property is also true for compressible flow. It has been shown by 
Cousteix and Houdeville in Ref. 36 that the behavior of solutions to the 3-D 
turbulent boundary layer equations is strongly effected by the hyperbolic 
character of the governing partial differential equations and by the nature of 
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the associated characteristic curves in the plane of the wing. As noted by 
them, when solved in the direct mode with the pressure distribution 
prescribed, the equations admit weak solutions in which the characteristics of 
the same family may Intersect each other on certain singular lines. The 
singular lines are analogous to shock waves across which the displacement 
thickness and other boundary layer properties are discontinuous. Such 
solutions have no physical significance and should not be identified with the 
separation lines observed in 3-D boundary layer separation. The jumps in 
displacement thickness leads to large (actually unbounded) values of the 
transpiration velocity which will act, through viscid/inviscid interaction 
mechanisms, to significantly alter the inviscid flow and streamline pattern 
near such singular curves. Since discontinuous solutions of this type cannot 
be self-consistent solutions of the coupled viscid/inviscid equations, it is 
clear that the interaction must eliminate the characteristic crossings and the 
associated discontinuities from the solution. One can speculate that the 
singular crossing curves are transformed through a strong local interaction, 
to an envelop of streamlines that could be identified with the locus of the 3- 
D separation. Cousteix and Houdevllle also proposed an inverse method for 
avoiding the discontinuous solutions of the direct problem. In their method, 
which was aimed at noninteracting type computations, they proposed specifying 
two Integral thicknesses and computing the pressure distribution as part of 
the solution of the boundary layer equations. 

Unfortunately, although Inverse methods may be useful for avoiding the 
shock line jumps In 3-D boundary layer solutions, they do not eliminate all 
problems faced in computing separated 3-D boundary layers. Indeed the major 
difficulty in computing such flows seems to be associated with the turning of 
the limiting streamlines and characteristics to a direction perpendicular to 
the marching direction and not with the formation of discontinuous weak 
solutions. In these cases the curved streamlines seem to form a "sonic" 
envelope across which the solution cannot be continued. This is so mainly 
because the flow upstream of such envelopes do not lie in the domain of 
dependence of the initial data line from which the solution was Initiated. 
Although the details are far from understood. It seems clear that such 
envelopes are likely to be relocated In some way to 3-D separation lines. The 
computation of 3-D boundary layers in these circumstances is extremely 
difficult because of the need to supply Initial data some place downstream of 
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the envelope and It Is not at all clear at present how this can be done. 
Because of the uncertainty in how to proceed in such cases, in our approach we 
employ a standard direct method to treat the viscid/inviscid coupling. 

If we assume there exists a solution of Eq (4.2.15) for the eigenvalues 
X = x i for i ■ 1,2,3, at any point (x,y n ) on the wing surface or in the wake, 
it follows that Eq (4.2.13) can be expressed in the form, 
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where , 
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M 3i * t-( b 2i “ a 2lV* b 12 " a 12 X 1* + (b ll “ a llV (b 22 " a 22 X i ^ /D i 

and 

D i “ ^ b 21 " a 2lV* b 32 “ a 32V " * b 31 " a 3lV* b 22 " a 22V 


for i * 1,2,3. 

Equations (4.2.16) are the compatibility equations of the hyperbolic 
equations in their "normal characteristic form". The solution of Eq (4.2.16) 
could be carried out most naturally with a "X type-scheme" (Ref. 45). For 
interior points, y 3 y n , the y derivatives of fi and y associated with 
each Xj factor are evaluated using values of etc. at (x, y n+1 ) and (x, y n ) 
when Xj < 0 and using values at (x,y n ) and (x, y n .^) when x^ > 0, where y n+1 > 
y n > y n _i* The x-derlvati ves of 6^, fl , and y then are computed using Eq 
(4.2.16) and new values of 6^, R , and y are obtained by an explicit 
marching scheme. There Is a fundamental problem that arises with any method 
of integration associated with the Implementation of the boundary 
conditions. The Integration domain lies between the wing root and tip 
stations. It Is seen that when Xj > 0 at y = y w1ng root or when x^ < 0 at y = 
y w ing tip* boundary conditions are required to supplement the compatabllity 
equations, Eq (4.2.16) at these stations. According to the rules in Kreiss 
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(Ref. 46), the total number of boundary conditions and compatibility equations 
should be equal to the order of the system of equations. The number of 
boundary conditions required at y * y^g root , * or exam P 1e > exactly the 
number of roots for which Aj >0. It is not clear at this stage how to 
properly determine these conditions. Our study indicates that computed 
results are very sensitive to the boundary conditions prescribed at the wing 
root and tip. For numerical convenience, in our method we employ zero 
spanwise gradient conditions. These conditions are non-physical and lead to 
errors in the region lying in the domain of dependence of the wing tip or body 
juncture points. Further study is clearly called for to determine correct 
procedures for setting these boundary conditions. 

4.2.3 Numerical Integration of Boundary Layer Equations 

In view of the difficulty in Imposing boundary conditions in the x- 
scheme, a less accurate but reliable numerical integration scheme is adopted 
for the solution of the system of Eq (4.2.11). The same procedure is also 
used to Integrate the laminar boundary layer equations. The method, first 
used by Smith, solves, first, at a constant x line, for the x derivatives of 
9^, R , y, and F. Then a two-level explicit integration scheme is used for 
integration in the x-direction. A local C.F.L condition Is imposed to 
determine the maximum integration step Ax. The y derivatives of the dependent 
variables e^, fl , y, and F are evaluated in such a way that the rules of 
domain of dependence and region of influence are not violated. Since the 
characteristic lines lie either on (as in the case of characteristics for F) 
or within an angle bounded by the external inviscld streamline (with an angle 
of inclination o with respect to the x-axis) and the limiting wall 
streamline (with an angle of inclination of a + e with respect to the x- 
axis), the use of these two directions to predict the directional bias of the 
disturbances always produces a conservative estimate of the step size to be 
used. We denote by (x, y n ) the location of the point where the y derivatives 
are to be evaluated and by (x, y n _i), (*» y n +i^ the two neighboring points 

where y n+1 > y n > y n _^. The y derivative |B. ^ (or, similarly for the 

3 6 

other derivatives - 5 -^ I , |J I v , |£- L w ) Is determined according 

<jy x,y n ®y *»y n *»J n 
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to the signs of the angles a and a + B* If both a and a + 0 are . 
positive, a backward difference of the form 


3fl I 
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form 


If both a and a + e are negative, a forward difference of the 


3fl , 

"57 'x,y n 


(R) 


x,y 


+ (H) 


n+1 
^n+1 " y n 


x,y r 


is used. And If a and a + e are of opposite sign, a central difference of 
the form 


3R , ,y n+l * ,Jf n-l 

^ x -*n ' y n+l - Vl 

is used. For the wing tip points, all the y derivatives are set equal to zero 
unless both a and a + 6 are positive and in this case a backward 
difference formula is used. For the wing root points, all the y derivatives 
are set equal to zero unless both a and a + e are negative and a forward 
difference formula is then used. 

The integration can be started with either laminar boundary layer 
equations or turbulent boundary layer equations at or near the leading edge 
lines of the wing. Because the Integration must follow the flow direction, 
the starting point must necessarily be downstream of the forward stagnation 
line. For calculations with a laminar boundary layer start, one must either 
Impose the location of transition to turbulent flow or use a natural transi- 
tion criterion to fix the transition position. In the present code we assign 
the transition point locations unless the laminar boundary layer solutions 
indicates separation, in which case we assign transition to the point of lam- 
inar separation. When transition occurs, the momentum thickness 8^ is 
assumed to be continuous across transition and the shape factor, fl , is 
assigned a value of 1.45 (an alternative to this is to assign a local jump in the 
value of R equal to Afl-|- rans jtj on as determined by experiments (Ref. 47)). 
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The 3-D boundary layer integration is extended into the wake in a manner 
similar to the 2-D calculation of Green et al. The skin friction is set 
equal to zero, and the dissipation length used in Green's formulation is set 
equal to half of that of the boundary layer value. Since Johnston's (Ref. 44) 
cross flow profile cannot admit a zero friction solution, only Mager's (Ref. 
43) crossflow profile can be used in the wake. In this case, y » e. 

Since the information of the characteristic lines is helpful in 
understanding the boundary layer flow behavior, the eigenvalues from the 
solution of Eq (4.2.15) are computed at each step of integration. For 
example, a very large value of X or appearance of imaginary roots of Eq 
(4.2.15) usually indicates the incipient breakdown of the computation. 

4.2.4 Special Considerations 

Since the solution procedure of the boundary layer equations as described 
in the previous sections is only an intermediate step of reaching the 
converged solution of the viscid-inviscld interaction, one does not always 
have a smooth invlscld pressure distribution for each cycle of the boundary 
layer computation. It is important in our Iterative method to establish a 
procedure for preventing a breakdown of the boundary layer computation from 
the appearance of an unphysical Intermediate inviscid pressure distribution. 
Our numerical experiments Indicate that a computational breakdown associated 
with large inviscid pressure gradients is often started by a rapid local 
growth of the shape factor FI . When this happens we set a maximum cutoff 
value of 2.4 for R and set the local y derivatives of the dependent 
variables equal to zero. The computed boundary layer solution of the shape 
factor R usually reaches a maximum at the trailing edge and drops off 
gradually going downstream into the wake. Occasionally, the computed value of 
R may drop to below 1 far downstream In the wake. These physically 
unrealistic values lead to numerical difficulties in the far field. To avoid 
such problems, we set a lower bound of R = 1.05 in the far field and 
continue the wake computation using a strip boundary layer approximation with 
all spanwise derivatives set to zero. 

Additional numerical difficulties may also arise because the streamwise 
momentum thickness, 8^, may occasionally drop below zero on the wing 
surface. The basic code would usually break down at such points because of 
appearance of a logarithm of in the computation of the skin friction from 
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the Ludwieg-Tillman formula. We have found, that in these cases the solution 
can usually be continued downstream using a strip boundary layer approximation 
with all spanwise derivatives set to zero. The program has been setup to do 
this if negative values of 9^ should arise at any stage of the computation. 

As mentioned in Subsection 4.2, we compute the metric coefficients 
hj^x.y), h 2 (x,y) and g(x,y) taking full account of the wing thickness and wake 
curvature. This more accurate treatment of the geometry turns out to be 
important near leading and trailing edges and acts to provide smoother, more 
accurate boundary layer solutions in these regions. To demonstrate this, we 
have computed two boundary layer solutions, one with our full metric 
coefficients expressions and the other with the metric coefficients for a 
corresponding flat wing. The computations are for the Lockheed Wing A at a 
Mach number M.. = .796, angle of attack a = 1.94 and reference Reynolds number 
of 5 million. The two solutions for the shape factor, R , are shown in Fig. 

9. It is seen that there are significant differences in the two solutions 
with the one using the full metric coefficients expressions showing milder 
gradients and a much lower peak value towards the trailing edge as compared 
with that of the flat wing solution. This behavior is generally true for the 
solution of all the other dependent variables. 

4.3 ITERATIVE SOLUTION TO VISCID-INVISCID INTERACTION ANALYSIS 

The inviscid and boundary layer solution procedures described in 
Subsections 4.1 and 4.2, must be modified to accomodate the matching 
conditions coupling the viscid and inviscid flows. In this subsection, we 
describe how the viscous coupling conditions are incorporated into the 
solution of the Inviscid equations and we also outline the interative 
procedure used to solve the coupled viscid/inviscid equations. 

4.3.1 Transpiration Coupling Conditions for Inviscid Flow Boundary 
Condi ti ons 

The boundary conditions for the finite volume formulation of the inviscid 
flow analysis need to be modified in order to take into account the boundary 
layer displacement and viscous wake effects. Transpiration boundary 
conditions are used In the present work. In 3-0 flow, the surface source 
velocity, m, is related to the displacement thickness, 6*, by the entrainment 
relation, Eq (4.2.7), which can be written in the form. 
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The displacement thickness, 6*, is defined by 


6 * = / 
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where pU Is the mass flux in the boundary layer in the external flow 
direction. The external quantities at the boundary layer edge, are Mach 
number (M), density (p e ), and the total velocity (u e ). The velocity 
components u, v, are components of u e along x, y, the boundary layer 
nonorthognal curvilinear coordinates. The quantities 1^, h 2 , are the metric 
coefficients and q is equal to 


q = [h*h| - 



1/2 


The displacement thickness Is determined as part of the boundary layer 
solution as described in Subsection 4.2. From Eq (4.3.1), therefore, the 
source, m, can be computed at each step of the boundary layer integration. 

The boundary conditions for the inviscid flow require the values of m at the 
midpoints of the inviscid mesh, which can be obtained by interpolation. Since 
the finite volume formulation is being used for the inviscid flow analysis, m 
is to be converted into M, its corresponding contravariant flux vector across 
the body surface. The scaling factor between m and M is quite complex and is 
not given here. It can be shown that for 2-D problem, M = m • ds , where ds 
is the arc length between the two nodal points of the Inviscid mesh. With the 
known distribution of the sources and its equivalent source jump in the wake 
(the difference of sources from top and bottom of the wake), the inviscid flow 
boundary conditions are modified as follows. 
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Referring to Fig. 10, the velocity potential, G, at the span station Z = 
Z K of the inviscid flow is to be computed. The value of G at Z = Z K+1 is 
lagged and computed from the previous iterate. The contravariant fluxes are 
to be computed at cell centers. On the boundary, A and E are the last points 
of the wake across the sheet, B and D are the trailing edge points and C is 
the leading edge point. On the wing surface (i.e., D - C - B) the reflection 
condition of the normal component of the contravariant fluxes, VM and VP, of 
the inviscid flow computation is modified due to the sources generated by the 
boundary layer. 


VM(I) + VP(I) =2-(p m + P p ) • (4.3.2) 

where p^ and pp are the densities computed at the two cells near the boundary 
and between the two space stations. The tangential components of the 
contravariant fluxes UM(I) = UP ( I ) , WM(I) = WP(I) etc remain the same as those 
of the noninteractive computation. On the wake (i.e., from E to D and from B 
to A) the flux conditions are not changed from the noninteractive computation 
because the image conditions have to be imposed, viz: 

UM(I) = - UP(M) 

VM( I) = - VP(M) (4.3.3) 

UM(M) = - UP( I ) 

VM(M) = - VP(I) etc 

where I and M are the respective image points across the wake sheet. The 
adding (or subtracting) of sources, however, is reflected in the mass flux 
balance equation (the continuity equation). Therefore, the residual formula 
evaluated at the nodes both on the wing and wake surfaces are modified to 
account for the mass injection from the boundary layer. Particular attention 
is paid to points D and B where there are source contributions both from the 
wing side and the wake side of the trailing edge. 

4.3.2 Make Curvature Effects 

In Subsection 4.1, we described the procedure for adjusting the mesh to 
follow the wake stream surface of the inviscid solution. In the inviscid 
computation, the wake surface was updated at each cycle of the inviscid 
Iteration. In the viscous computation It was found to be more efficient to 
update the wake position before each boundary layer computation rather than at 
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each inviscld cycle. In the viscous computation the wake position is defined 
in terms of a wake angle, f? , 



Z=Z K 


7 

0 



where X, Y are the Cartesian wake coordinates, & is the wake line angle with 
respect to X-axis at constant Z = Z K , 0 and 7 are the inviscid velocity 
components, averaged between the top and bottom of the wake surface, along X 
and Y-axes, respectively. 

The wake curvature effect terms are imposed in the following manner. 
Referring to Fig. 10, let the row of computation nodes at the wake, above the 
wake and below the wake at a spanwise station, Z = Z^, be denoted by ( I ,KY ,K) , 
(K,KY-1,K) , (I,KY+1,K), respectively, and let (M,KY,K) etc. be the image nodes 
across the wake. The Kutta condition and the wake curvature condition can be 
Imposed according to the following formulae, for the surface wake nodes, 

G(M,KY,K) = G(I,KY,K) + CIRC(K) + T (4.3.4) 

and for the image wake nodes, 

G(M,KY+1 ,K) = G(I,KY-1,K) + CIRC(K) + r (4.3.5a) 


G(I,KY+1,K) = G(M,KY-1,K) - CIRC(K) - r (4.3.5b) 


The function G is the reduced velocity potential, CIRC(K) is the trailing edge 
velocity potential jump determined by the Kutta condition and r, the 
circulation, is equal to the jump in velocity potential across the wake. The 
circulatory function, r, is determined by the matching conditions coupling the 
viscid and inviscid flows. Within standard boundary layer theory it can be 
written in the form. 


r(s) | 

Z * Z K 
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te 


(4.3.6) 
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where 0 is the average surface value of the outer inviscid velocity in the 
chordwise direction, and are the respective sums of the upper and lower 
displacement and momentum thicknesses, respectively. Our method incorporates 
corrections to the standard coupling conditions to account for strong 
interaction effect at trailing edges. These corrections modify the standard 
expression for r given above, as described in the following section. 

4.3.3 Trailing Edge Corrections 

In our method we incorporate strong interaction effects at trailing edges 
into our basic vlscid-inviscid coupling procedures. We follow closely the 
procedures developed for the 2-D airfoil problem as described in Ref. 6. In 
that work a local asymptotic solution is developed which describes the strong 
interaction effect at trailing edges. Through the use of this solution, 
modifications to the classical coupling conditions are developed to take into 
account normal pressure gradient effects across the boundary layer in the 
trailing edge region. The theory of Ref. 6 was for a strictly 2-D flow. In 
the present work, we adapt the 2-D corrections of Ref. 6 to the wing problem 
considered in this study using a quasi -2-D approximation. In this approach, 
we simply apply the 2-D corrections to the full 3-D form of the classical 
coupling conditions assuming the flow Is 2-D in streamwise planes. The use of 
a 2-D strip approximation is justified because the flow gradients normal to 
the trailing edge are asymptotically larger compared to gradients in the 
spanwise direction and the local analytic solution will be 2-D under these 
circumstances except the near wing tip and body juncture stations. To be 
consistent for swept trailing edges, the quasi -2-D theory should be applied in 
a direction normal to the trailing edge. However, for simplicity we apply the 
theory in the streamwise direction. Since the corrections are relatively 
small we believe the simple theory will be adequate except, perhaps, for very 
highly swept trailing edges. 

In the modified theory, corrections are applied to the standard coupling 
conditions to account for normal pressure gradient effects. The corrections 
are applied to 1) the source velocity on the wing surface, 2) the jump in 
source velocity in the wake and 3) the jump in pressure or velocity potential 
(i.e., D across the wake. In addition, the surface pressure computed from 
the inviscid solution is corrected to account for the pressure drop across the 
boundary layer in the trailing edge region. The corrected boundary 
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conditions, following from Ref. 6, are written in the form: 

Source velocity: 

m(s,Z K ) = m o (s,Z K ) [/(s-s^ K( s .Z^) ] wing section (4.3.7a) 


[m(s,Z K )] = Cm + ( s te ,Z K ) R + (s,Z K ) -m"(s te ,Z K )R"(s,Z K )] + m w (s,Z K ) 


wake 


(4.3.7b) 


Wake Circulation: 


r ( s ) *Z=Z = “ 1 / 0 («* + 0 n ) A s-s te ) J w ( s)ds] z=z (4.3.7c) 

K £ K 

0 te 

where m 0 (s, Z^) and m^s, Z^) are the classical expressions for the source 
velocity on the wing at section Z = Z K and the jump in source velocity in the 
wake, m + (s te , Z K ) and m" (s te , Z^ ) are the values of the corrected source 
velocities on the upper and lower surfaces of the wing section at the trailing 
edge, s is the arc length along the wing section and wake. The functions K(s, 
Z K ), R(s,Z k ) and J w (s) are given in terms of universal functions from the 
local trailing edge solution in Ref. 6. The corrected pressure distribution 
in the wake is given by an expression of the form 


P e (s ’ Z K ) = P e (s, V " TTa t p ] 


(4.3.8) 


where P g (s,Z K ) is the corrected pressure distribution on the wake, 

P*(s,Z.,) is the inviscid solution for the pressure on the upper surface of the 
wake and [ P ] is the difference in pressure across the wake from the Inviscid 
solution. The jump in pressure in the inviscid solution is caused by the wake 
boundary condition Imposing the jump in velocity potential, r, across the 
wake. The quantity A is also given in terms of universal functions in Ref. 

6. The expression for the pressure distribution on the wing section surface 
Is written in the form, 

P e (s, Z K ) . P e (s, Z K > - [ P ] te Id. Z K ) - P e Q| p 0 d, Z K > (4.3.9) 
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where ^(s, Z^) and P e (s, Z^) are the corrected and inviscid surface pressure 
distributions respectively, [ P ] is the jump in inviscid pressure across 
the wake evaluated at the trailing edge. The quantity I ( s, Z^) is a 
universal function given from the local trailing edge solution and p 0 (s, Z^) 
is the pressure correction determined from the classical boundary layer theory 
as described in Ref. 6. 

4.3.4 Global Iteration Strategy for the Interaction Solution 

The solution procedure for the global iteration is carried out in the 
following steps, 

I. Obtain the inviscid 3-D solution using the AFZ scheme, compute the 
surface velocities and surface flow angles. 

II. Compute the corrected surface velocities and surface flow angles as 
indicated in the previous subsection. 

III. Obtain the 3-D boundary layer and wake solution using the modified 
surface velocities and flow angles from the previous step. 

IV. Compute the coupling condition as follows: 

1. Compute the surface source velocity distribution and its 
contravariant equivalent from the boundary layer solution. 

2. Compute the floating wake surface coordinates using the Inviscid 
velocities, and obtain the equivalent shear function used in the 
inviscid parabolic coordinate mapping. 

3. Compute the wake jump condition from the updated wake 
coordinates. 

V. Go to step I. 

There are two sets of relaxation factors assigned to the calculations, 
the P^, 1 = 1,2,3, for the inviscid solutions and the relaxation factors V^, i 
■ 1,2, 3, 4, for the boundary layer source flux on the wing, the viscous wake 
source flux, the floating wake coordinates and the velocity potential jump for 
the wake curvature effect. The chosen values of P^ are based solely upon the 
best convergence behavior of the inviscid calculation. The values of range 
from .5 to 1. Since the inviscid AFZ scheme requires four sweeps per one 
outer cycle of iteration, the boundary layer calculations are updated after 
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multiples of four sweeps of the inviscid iteration. Eight sweeps of inviscid 
iterates are chosen between boundary layer calculations for our computations. 

Our experience indicates that the inviscid part of the program is quite 
robust and convergence has been achieved for quite a broad range of test cases 
including wings of practical interest. The major requirement for a successful 
convergence calculations for the viscous wing program during a global 
iteration is the avoidence of boundary layer computation breakdown. For a 
modern transonic wing with large trailing edge camber and wing twist, large 
inviscid pressure gradients can cause a breakdown of the boundary layer 
computation at the outset. Since viscous effects generally smooth out these 
large pressure gradients, such breakdowns in the initial stages of the 
computation do not necessarily imply that a converged solution cannot be 
obtained. Some practical techniques can be used to achieve convergence. For 
example, if the computational problem is the result of too strong an inviscid 
shock locally, the Iteration can be initiated at a lower freestream Mach 
number to achieve an intermediate convergence, then the Mach number is 
progressively increased to the desired value. On the other hand, if the 
difficulty in obtaining a boundary layer solution is due to the complexity of 
the camber variation and wing twist, it is desirable to start the iteraction 
at a low angle of attack which is then increased to the desired value. These 
techniques shall be implemented in future versions of the code. 

With a given wing geometry, a given case requires the specification of 
the free stream Mach number M,,, the angle of attack, a, and the reference 
Reynolds number Re re f. Numerical computation is performed with an assigned 
mesh and sets of relaxation parameters Pj and V i . In addition, the boundary 
layer transition location needs to be assigned. An improperly assigned 
transition location can lead to erroneous results, as shall be discussed in 
detail in the next section. A typical converged interaction calculation for a 
supercritical flow case takes about 10 to 20 boundary layer cycles of 
computation dependent on the relative difficulty of the case. Normally it 
takes about eight cycles of boundary layer computation to reduce the maximum 
Inviscid residual to less than 10" 5 . 

The convergence history of a typical case run for the Lockheed Wing A is 
shown in Fig. 11. The free stream Mach number used was M a = .82, the angle of 
attack a - 1.5°, reference Reynolds number Re re f = 18 x 10 6 , and grid 
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distribution (160 x 16 x 32). The wing root circulation value and the total 
number of supersonic points are plotted as a function of iteration numbers. 
The CPU time required for the total thirty cycles of boundary layer 
calculations is about eight minutes using a Cray-lS computer. 


50 





















5. RESULTS 


5.1 NUMERICAL RESULTS - INVISCID FLOW 

Results for the ONERA M6 wing of Fig. 12 are illustrated in Fig. 13 
through 22. For M* = .923, a = 0, the computed C p distribution is given in 
Fig. 13. It can be seen that there is a large supersonic zone with shocks 
extending over the length of the wing. A comparison with experiment near mid- 
span is given in Fig. 14. From Fig. 15 it can be seen that the solution is 
well converged after 32 iterations. Figure 16 compares the convergence 
history of the average residual for two variations of the present AFZ method 
with a relaxation (SLOR) method. The development of the supersonic zone for 
the three methods is shown in Fig. 17. Figure 16 seems to indicate that SLOR 
exhibits a faster convergence at the beginning of the iteration. This effect 
is even more pronounced in the lifting case to be considered later (Fig. 

21). The convergence of SLOR slows down significantly after the first 50 
iterations while the AF schemes maintain a steady and rapid convergence 
rate. The advantages of the AF scheme shows up even more clearly when 
considering the convergence of the global aspects of the solution such as the 
number of supersonic points (Fig. 17) and the lift (Fig. 22). 

Figures 18 through 22 give results for the same wing (Fig. 12) at M,, = 

.84 and o= 3.06°. Figure 18 shows the upper and lower surface pressures. It 
can be seen that the two shocks (one near the leading edge, the other further 
back) come together as the wing tip is approached. Figure 19 shows a 
comparison with experimental surface pressures, at a spanwise section near 
mid-span. Figure 20 shows the convergence in surface pressure. Figure 21 
shows the convergence history in terms of average residual. Again it seems 
that SLOR has a faster Initial convergence rate, but considering Fig. 22 
(development of lift) the advantage of the AF schemes is obvious. 

The computation of the invlscid transonic flow about transport type wings 
like the ONERA M6 is relatively straightforward. This is because these wings 
are characterized by high aspect ratios, low sweep, small twist and similar 
airfoil sections. The 3-D effects in the flow field about these types of 
wings are minimal. The results of Fig. 23 and 24 are intended to demonstrate 
the capability of the inviscid code for computing the flow about highly 3-D 
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wings. Figure 23a shows a 3-0 view of a low aspect ratio, highly swept and 
twisted wing. The flow field about this wing was computed at a free stream 
Mach number, M,,, = .9 and an angle of attack, a = 8°. The surface pressure 
distribution at three sections are shown in Fig. 23b, 23c, and 23d. Figure 
23b shows the root airfoil section and the corresponding surface pressure. 

The section lift at this station is small because of the lower relative 
incidence. Figure 23c shows computed results at the mid-span section. The 
section load here is higher than that at the root although the section is 
twisted down relative to the root. This effect is due to the upwash from the 
root section caused by the large sweep of the wing. This effect Is even more 
pronounced at the tip section (Fig. 23d) where the large relative incidence 
(due to upwash) causes a strong shock at the leading edge. 

Figure 24a shows the planform and surface Isobars for another wing. This 
wing was also run at a Mach number of M,* = .9 and o= 8°. The isobars clearly 
show a two shock pattern before the crank in the trailing edge. The leading 
edge and trailing edge shock waves Interact near the wing tip. Figure 24b 
shows the airfoil section and surface pressure at the root of the wing. 

Figure 24c shows the mid-span section where the two shocks can be seen 
clearly. Figure 24d shows the tip section after the two shocks have 
interacted. The leading edge camber of the outboard sections of this wing is 
shown in Fig. 24d. The effect of leading edge camber, in eliminating strong 
leading edge shocks, is demonstrated dramatically by comparing Fig. 24d and 
23d. Both these wings are highly 3-D and heavily loaded at their tips. These 
conditions usually require adjustments in the parameters of the AFZ inviscid 
code. For these more difficult cases the iteration procedure should be 
underrelaxed somewhat, particularly for fine grids. The fact that converged 
solutions could be achieved in these very difficult cases indicates the 
soundness of our basic AFZ scheme. 

A comparison of convergence rates using the two far field boundary 
conditions mentioned in Subsection 3.1 Is shown in Fig. 25. In the figures 
the solid line was computed using the expansion of the Prandtl -G1 auert 
equation (Eq 3.1.8) the dashed line was computed using the original condition 
4> = 0 upstream and = 0 far downstream. Figure 25a shows how convergence is 
improved using the new far field condition. Figure 25b shows that the lift is 
increased significantly with the new condition. The computed lift using our 
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Prandtl -G1 auert far field boundary condition Eq (3.1.8) seems less sensitive 
to the location of the far field boundary (i.e., less sensitive to the 
parameter BOUND). The results of Fig. 25 were computed with a BOUND of .95 
which places the outer boundary about 6 cords from the wing. 

Figure 26 shows a comparison of computed results using the first and 
second order artificial viscosities. Figure 26a shows these convergence 
rates. The figure shows that the convergence is slowed using the second order 
scheme. Figures 26b-26j shows sectional surface pressure distributions from 
the root to tip. The figures show that the pressure distributions are all 
sensitive to the choice of the artlfical viscosity formula. The shock in the 
second-order results is sharper (for example, see Fig. 26g) and the suction 
peak near the leading edge is significantly reduced. The lift is increased 
from = .635, first order, to = .665, second order. 

5.2 VISCOUS INTERACTION RESULTS 

In this subsection, we present results computed from the viscous wing 
program. The computer program developed for the present program is designated 
GRUMWING. The present version of the code is a pilot code which we expect to 
further develop into a production code. The inviscid results presented 
earlier In this Section were computed with the inviscid option of the GRUMWING 
program. Details of this program are given In the user's manual Included as 
the second volume of the report. 

The results presented in this section were obtained for the Lockheed 
Wings A and B described in the report of Hinson and Burdges (Ref. 48). Of 
particular significance in their work is the fact that special attention had 
been paid to the effects of the fuselage on wing pressure data so that 
meaningful comparisons can be made for an analysis method with wing alone. 
Shown in Fig. 27 and Fig. 28 are the wing planform and wing airfoil sections 
for the Wing A and the Wing B, respectively. Wing A Is a transonic wing of 
transport type. It has an aspect ratio of 8, and quarter chord sweep of 
25°. Wing B is a fighter type wing designed for transonic cruise. It has an 
aspect ratio of 3.8, and quarter chord sweep of 30°. The nominal test 
Reynolds number for the experiments based on the mean aerodynamic chord was Re 
= 6 x 10 6 for Wing A and Re = 10 x 10 6 for Wing B. 
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5.2.1 Boundary Layer Transition Assignment 

It has been known that Inappropriately selected transition points can 
lead to unrealistic results even when the computation proceeds without any 
apparent difficulty. In this section we show results of computations with 
varying transition point locations that give some Indications of the 
sensitivity of the solution to the choice of transition point location. 

In order to demonstrate the effects of the transition point locations on 
the solution, computations were carried out for a variety of transition 
locations for the Lockheed Wing A at = .82, a = 1°, Re = 5 x 10® using a 

(160 x 16 x 32) grid. Shown in Table 1 are the transition results for runs 
using assigned transition locations of 0, 2.5, 5, 7.5, and 10% chord, on the 
upper and lower surfaces. For the last case laminar separation occurs on the 
upper surface at 8.2% chord. For this case, the transition specification is 
overridden and transition is set at the laminar separation point. For all 
cases, turbulent separation occurred on the lower surface which Is physically 
unrealistic and is due to the assignment of the transition point too far 
forward on the lower surface. To show this, four runs were made for the same 
Mach number, angle of attack, reference Reynolds number and transition point 
location (10% chord) on the upper surface. Calculations were run with 
transition locations on the lower surface of 20, 40, 60, and 80% chord. 

Results are shown in Table 2. It Is seen that all four cases end up with 
transition on the upper surface set at 8.2% chord (laminar boundary layer 
separation) Induced by the large suction peak in the pressure distribution. 

The turbulent boundary layer on the upper surface remains attached. The 
solution on the bottom surface separates at 88.2% when the transition point is 
set at 20% chord. For transition point at 40% the flow remains fully 
attached. In the remaining two cases of assigned transition at 60 and 80% 
chord, the transition occurs naturally at 47.8% chord due to laminar boundary 
layer separation resulting In identical solutions with no turbulent 
separation. 

From the above discussion, it is clear that the solutions can be 
sensitive to the transition point location, particularly when the transition 
point on the lower surface is set far forward. For a typical transonic wing 
geometry, the transition point on the upper surface should be set just down 
stream of the suction peak. For the lower surface, if the computation is made 
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for wind tunnel data comparisons, the separation point should be set at the 
location where the boundary layer Is tripped, otherwise, the transition 
location should be set at far enough downstream of the leading edge t6 Insure 
natural transition. 

5.2.2 Reynolds Number Dependence of Viscous Wing Solution 

Computations were made to study the Reynolds number dependence of the 
theoretical solutions. The computations are for the wing A at M,,, = 0.82 and a 
= 1.5°. Reference Reynolds numbers for the seven cases were, 5M (5 million), 
8M, 10M, 12M, 14M, and 16M, respectively. The runs above Re = 10M required 
smaller integration steps in the boundary layer solution in order to avoid 
unrealistic results associated with the momentum thickness becoming 
negative. The results for the lift variation with Reynolds number are shown 
In Fig. 29. There Is a strong Reynolds number effect for Re < 8 x 10®. The 
lift levels off quite rapidly beyond 10M indicating that the results of free 
flight chord Reynolds number in the 100M - 200M range can be extrapolated from 
the present results. The inviscld, large Reynolds number limit for the lift 
coefficient in this case is C L * .637 as determined from an inviscld 
calcualtion. These results show that very large Reynolds numbers are required 
to approach the Inviscid limit and that viscous effects will be very important 
even at flight Reynolds numbers. 

5.2.3 Solution for Pressure & Section Lift Distributions & Comparison with 
Experiment 

In this section, we compare solutions for the pressure distribution, 
lift, and drag with experimental data for Lockheed Wings A and B. The 
pressure distributions given here are the composite pressure as given by Eq 
(4.3.8) and (4.3.9). The freestream Mach number, M^,, and chord Reynolds 
numbers were matched with those of Lockheed test cases. The boundary layers 
were tripped at 5% chord from the leading edge. Since the effective angles of 
attack of the wind tunnel test are subjected to major uncertainties due to 
wall Interference, we carried the computations out at an angle of attack 
chosen to closely match experimental values of total lift. 

For the Wing A run, the free stream Mach number M*, was .82, the angle of 
attack a was set to 2°, the reference Reynolds number Re re f was six million 
and a grid of 160 x 16 x 32 was used. In addition to a 5% chord tripped 
boundary layer run, a natural transition run was made for comparison. 
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Furthermore, an inviscid run was also made for the same M,, and a. Nothing 
special was required to achieve a converged solution. The maximum inviscid 
residual was down to less than 10"^ after about 10 cycles of boundary layer 
computations or ninety inviscid fine grid sweeps. At this point the Inviscid 

pressure was accurate to within 10“^ and CPU time was about three minutes on a 

Cray-lS computer. Of this about 2 minutes of CPU time was for the computation 
of the inviscid solution. The convergence history of the total number of 
supersonic points and the wing root section circulation are shown in Fig. 

30. The position of the floating wake was the last quantity to converge. 
Figure 31 shows the convergence history of the last point of the wake at the 
wing root section as a function of boundary layer iteration. The converged 
wake at this section is shown in Fig. 32. 

The results show that for the natural transition case, the transition 
occurred at 8.1% chord on top of the wing and 50% chord on the bottom 

surface. The results for the forces are shown In Table 3. It is seen that 

the wing lift coefficients, C L , and drag coefficients, C D , for the two 
transition calculations are relatively close to each other. The two transi- 
tion calculations gave different boundary layer solutions and hence different 
friction drags. This will be discussed in the next section. The C L for the 
viscous solution (tripped boundary layer) was .536 as compared with the 
inviscid value of .749. The experimental value of C L was .53 at a nominal 
angle of attack of 2.9° (compared to the a = 2° used in the comutation). 

Results for section pressure distributions (for the tripped boundary 
layer case) are shown in Fig. 33a - 33e for n = .15, .30, .50, .70, and .95, 
respectively, where n is the sectional distance from the wing root normalized 
with the half span. Shown also are the Inviscid pressure distributions at the 
same free stream Mach number (M,,, = .82) and angle of attack (a = 2°) as well 
as the experimental results. It is seen that the viscous effect at this 
Reynolds number (6M) is quite strong. Shock location is shifted upstream and 
shock strength becomes much weaker resulting in a much lower lift. The 
viscous effect on the bottom wing surface is relatively small. The reasons 
for the poor agreement of the pressure distribution with the data on the lower 
surface near the trailing edge are not understood at this time. The 
comparisons for span load distributions for the same case are shown in Fig. 

34. We note again the vast differences between the viscous and inviscid 
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solution for the section lift distribution. The isobar plots for the 
converged viscous solution are shown in Fig. 35. Pressure distributions on 
the wing and wake for the top and the bottom surfaces are displayed. Because 
of the viscous effects, the shock patterns differ greatly from those in the 
inviscid solution. The lambda type shock pattern is clearly exhibited. The 
bottom surface shows relatively smooth flow patterns with no shock waves. The 
pattern is "2-D" like on the lower surface except near the wing tip region. 

For the Wing B runs, the free stream Mach number was M,,, = .9, the 
Reynolds number Re re f was 10 million and the boundary layers were tripped at 
5% chord to match the experimental values. The computation grid was (128 x 16 
x 32). Two angles of attack a = 3° and o = 3.5° were run in order to compare 
with the nominal experimental value of a= 3.9°. The results for section lift 
distributions at n = .2, .4, .6, .8, and .95 and wing loadings are shown in 
Table 4. Because of higher Reynolds number, it is seen that the lift results 
for the viscous flow calculation differ by only about 10% from the inviscid 
values. The experimental value of C L was .49 at a nominal angle of attack of 
a * 3.9°. Sectional pressure distributions for both the viscous and inviscid 
calculations at a = 3.5° are shown in Fig. 36a - 36e for n = .2, .4, .6, .8, 
and .95, respectively. The span load distribution for the same case is shown 
in Fig. 37. The overall viscous effects are small as compared to the Wing A 
computations which is consistent with the difference in Reynolds numbers in 
these two cases. 

5.2.4 Boundary Layer & Wake Solutions 

In this section, we present solutions for some of the boundary layer 
quantities from the results of the converged viscid-inviscid interaction 
solution. For Wing A calculations at M,* = .92, a = 2° and Re re j: = 6 x 10 6 , 
two solutions are presented, one for tripped transition (5% chord) and one for 
natural transition. The results for the displacement thickness <5*, the 
momentum thickness 0^, the skin friction and the shape factor fl at the 
mid span station are given in Fig. 38 - 41, respectively. In the natural 
transition case, the transition point on the upper surface occurred at 8.1% 
chord. The two solutions on the upper surface are virtually identical. On 
the bottom surface, natural transition occurred at 50% chord compared to the 
5% location of the fixed transition case. The solutions on the lower surfaces 
differ significantly. The displacement and momentum thicknesses are much 
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smaller in the natural transition runs as shown in Fig. 39. The comparison of 
the skin friction results shown in Fig. 40, indicates that natural transition 
will result in a smaller friction drag due to a larger area of the wing 
surface covered by laminar flow. 

The isocline results for the tripped transition run are shown in Fig. 42a 
- 42f. Both the top and the bottom wing and wake surface plots are shown. 

The isoclines for the surface velocity are shown in Fig. 42a. This velocity 
is the converged interaction "inviscld solution" with trailing edge strong 
interaction corrections. It is seen that the contours are quite similar to 
those of the isobar results shown in Fig. 35. On the top wing surface, the 
lambda shock structure is restricted to about 40% of the half span region from 
the wing root. In this region, the boundary layer exhibits more 3-D 
effects. On the bottom wing surface, the spanwise variations of the solution 
is relatively small. The wing tip effects are apparent as shown in the 
results of surface velocity angle a in Fig. 42b. Relatively large gradients 
are seen near the tip and towards the trailing edge and into the wake. It is 
also seen that the contours of the surface velocity angle in the wake for the 
top and the bottom surfaces are quite similar indicating that the viscous wake 
in this case does not have a very complicated 3-D structure. The boundary 
layer solutions for the displacement thickness 5*, the momentum thickness e^, 
the shape factor fl and the limiting wall streamline angle 3 are shown in 
Fig. 42c - 42f, respectively. The displacement thickness 6* has a rather 
curious wake distribution on the top surface where a moderate gradl-ent in the 
spanwise direction is shown near the wing root region. This reflects the fact 
that more realistic wing root boundary conditions are probably required. This 
aspect is beyond the scope of the present study. From our experience, the 
shape factor fl is the most sensitive part of the boundary layer 
calculation. At early stages of the interation, there is relatively large 
spanwise gradient for fl near the intersection of the two shocks. Boundary 
layer computations can break down if the shocks are too strong. At a later 
stage, the spanwise gradient gradually decreases as the viscous effect sets in 
and shocks become weaker. The converged results for fl is shown in Fig. 

42e. The distribution of the limiting wall streamline deviation angle 3 from 
the external inviscid streamline is shown in Fig. 42f. 

The variation of 3 is relatively mild except near the outboard and trailing 
edge region on the bottom surface, where most of the spanwise flow occurs. 


58 



6. CONCLUSIONS 


The principle conclusions of this report are: 

o The approximate factorization scheme (AFZ) used in the inviscid 

calculation presented in this report is as efficient as any state-of- 
the-art scheme available today. Inclusion of the Prandtl-Glauert far 
field boundary condition made the scheme faster and more reliable. 
Extending the method to second-order accuracy in the supersonic zone 
affected the solution significantly in critical areas (l.e., leading 
edge and wing tip). The assumption of irrotational flow does 
introduce some inaccuracies, particularly for strong shocks. The use 
of the Euler equations may be an alternative but their solution would 
require more computing time and memory. In addition, the Euler 
equations may Introduce additional computational difficulties. 

o The 3-D boundary layer equations resulting from the integral method 
formulation are very complex. There Is room for improvement in the 
solution method of these equations. The lack of a complete set of 
physical boundary conditions at the tip and root sections forces us 
to use arbitrary zero gradient conditions at these stations and this 
Introduces errors on the wing in the region of influence of these 
points. For high aspect transport wings, errors introduced In these 
regions are small and the consequences are not too important for 
engineering purposes. However, this will not be the case for highly 
swept, low aspect ratio fighter type wings. 

o The results Indicate that the wing loadings are strongly controlled 
by the inviscid solutions which are affected by the grid resolution 
and distribution. 

o The "direct mode" boundary layer analysis as employed in the present 
study has a supri singly large region of validity. While the "Inverse 
mode" boundary layer analysis has proved very successful in 2-Q 
studies with small separated region, the method has yet to be 
explored in 3-D boundary layers. Three-dimensional boundary layer 
separation is Inherently more complex and very likely cannot be 
adequately treated by extension of the 2-D Inverse methods. 
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o The present investigation has demonstrated that the "zonal" approach 
is a practical procedure for computing viscous flows over realistic 
wings at transonic speeds. Results are qualitatively good in many 
cases although many areas of Improvement are required to achieve 
accuracies comparable to our 2-D viscous airfoil codes, GRUMFOIL. 
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Fig. 8 Coordinate System for 3-D Boundary Layers 
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Fig. 17 Development of Supersonic Zone, Nonlifting Case 
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Fig. 19 Comparison with Experiment, Lifting Case 




















NO. ITERATIONS 


Fig. 21 Convergence Comparison: Average Residual, 

Lifting Case 
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Fig. 23 Sample Wing7¥ = 6.9, a * 8° (Sheet 1 of” 2) 
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c) Midspan Section 



d) Tip Section 

Fig. 23 Sample Wing, * 0.9, a - 8° (Sheet 2 of 2) 
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a) Planform & Isobars 



b) Root Section 

Fig. 24 Sample Wing 2, = 0.9, a = 8° (Sheet 1 of 2} 
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Fig. 25 Comparison of Far Field Boundary Conditions 
Lockheed Wing A, = 0.82, o = 1.5° 
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d) Z = 4.5 

Fig. 26 Comparison of 1st and 2nd Order Artificial 

Viscosity Lockheed Wing A, M = 0.82, a = 1.5° 
(Sheet 2 of 5) 
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f) Z = 9.0 


Fig. 26 Comparison of 1st and 2nd Order Artificial 

Viscosity Lockheed Wing A, M = 0.82, a = 1.5° 
(Sheet 3 of 5) 
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h) Z = 13.5 

Fig. 26 Comparison of 1st and 2nd Order Artificial 

Viscosity Lockheed Wing A, M = 0.82, a = 1.5° 
(Sheet 4 of 5) 
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j) Z = 18.0 


Fig. 26 Comparison of 1st and 2nd Order Artificial 

Viscosity Lockheed Wing A, M m = 0.82, a = 1.5° 
(Sheet 5 of 5) 


93 












PRESSURE 

DATA 

SECTIONS 



WING A 


TJ*0 



WING REFERENCE PLANE 


i? - 0.5 


17 * 1.0 


WING A AIRFOIL SECTIONS IN RIGGED POSITION 
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Fig. 28 Lockheed Wing B 
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Fig. 33 Comparisons of Sectional Pressure Distributions 
for Lockheed Wing A (Sheet 1 of 3) 
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Fig. 33 Comparisons of Sectional Pressure Distributions 
for Lockheed Wing A (Sheet 2 of 3) 
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99 





Fig. 34 Comparisons of Span Load Distribution for 
Lockheed Wing A 
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Fig. 36 Comparisons of Sectional Pressure Distributions 
for Lockheed Wing B (Sheet 2 of 3) 
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Fig. 37 Comparisons of Span Load Distribution for 
Lockheed Wing B 



Fig. 38 Displacement Thickness Distribution at 
Mid-Half-Span Station 




























Fig. 39 Momentum Thickness Distribution at Mid-Half-Span 
Station 
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Fig. 40 Lower Wing Surface Skin Friction Distribution 
at Mid-Half-Span Section 
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c) Displacement Thickness 
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Fig. 42 Converged Isoclines for Viscous Wing Solution 
(Sheet 2 of 3) 
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TABLE 1 TRANSITION RUN RESULTS USING LOCKHEED WING A 
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TABLE 2 TRANSITION RUN RESULTS USING LOCKHEED WING A 
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TABLE 3 COMPUTED FORCES RESULTS FOR LOCKHEED WING A 
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TABLE 4 COMPUTED FORCES RESULTS FOR LOCKEEO WING B 


Section Lift C L 

n 

Viscous 

Solution 

Invi scid 

Solution 

Experiment 

mam 

3.5° 

3° 

3.5° 

3.9° 

0.2 

0.4516 

0.5024 

0.4966 

0.5539 

0.48 (n =-22) 

.4 

0.4738 

0.5301 

0.5252 

0.5882 

0.52 

.6 

0.4753 

0.5351 

0.5339 
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0.55 

.8 

0.4458 

0.5037 

0.5116 

0.5746 

0.51 

.95 

0.3588 

0.4056 

0.4229 

0.4771 

0.38 


WING LOADING 


c. 


c n 

Pressure Drag 

Viscous 

Invlscid 

Viscous 

Invi scid 

a = 3° 0.4483 

0.5006 

0.0190 

0.0217 

a = 3.5° 0.5020 

0.5605 

0.0242 

0.0280 

a = 3 . 9 ( EXPERIMENT ) 0.49 



— 

= 0.9 Re ref = 10 M 

Ak(l) = 

Ale (2) = 0.05 

Grid = 128 x 16 x 32 
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one scalar version to be used on an IBM-3081 and two vectorized versions on 
Cray-1 and Cyber-205 computers. 
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